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;
using Microsoft.VisualBasic;

namespace vcsSWG
{
    struct DRAWPARA
    {
        public int kpt;            //拡大率（画面表示：1，ファイル保存：2）
        public int HSIZE;          //画像横寸法（px）
        public int DHL;            //画像左余白寸法：縦軸描画用（px）
        public int DHR;            //画像右余白寸法（px）
        public int VSIZE;          //画像縦寸法（px）
        public int DVL;            //画像下余白寸法：横軸描画用（px）
        public int DVU;            //画像上余白寸法（px）
        public string sxjiku;      //ｘ軸名（単位含む）
        public double xmin;        //ｘ軸最小値
        public double xmax;        //ｘ軸最大値
        public double dx;          //ｘ軸増分
        public string slogx;       //ｘ軸の種類（L：対数，N：普通）
        public string syjiku;      //ｙ軸名（単位含む）
        public double ymin;        //ｙ軸最小値
        public double ymax;        //ｙ軸最大値
        public double dy;          //ｙ軸増分
        public string slogy;       //ｙ軸の種類（L：対数，N：普通）
        public int kxi0;           //グラフ左端座標
        public int kxf0;           //グラフ右端座標
        public int kyi0;           //グラフ上端座標
        public int kyf0;           //グラフ下端座標
    }

    public partial class Form1 : Form
    {
        const double pi=Math.PI;
            
        public Form1()
        {
            InitializeComponent();
        }

        //---------------------------------------------------------------------------------
        private void Form1_Load(object sender, EventArgs e)
        {
            toolStripLabel1.Text = "繰返回数";
            toolStripLabel2.Text = "誤差（目標SP定義範囲）";
            toolStripLabel3.Text = "誤差（全体）";
        }
        //---------------------------------------------------------------------------------
        private void toolStripButton1_Click(object sender, EventArgs e)
        {
            toolStripLabel1.Text = "繰返回数";
            toolStripLabel2.Text = "誤差（目標SP定義範囲）";
            toolStripLabel3.Text = "誤差（全体）";
            main_part();
            MessageBox.Show("計算完了", "完了通知");
        }
        //---------------------------------------------------------------------------------
        private void main_part()
        {
            double dt;       //時刻歴波形サンプリングピッチ
            int ndata;       //時刻歴波形データ数
            int nn;          //FFTにかけるデータ数（２の累乗）
            int nfold;       //折り曲げ振動数を与える番号に１を加算した数値
            double h=0.05;   //減衰定数
            double ddymax;   //最大加速度
            double dymax;    //最大速度
            double ymax;     //最大変位
            double eps1=0.0; //目標スペクトル定義周期範囲での誤差
            double eps2=0.0; //全誤差
            double rk;       //目標スペクトルと模擬地震波応答スペクトルの加速度比
            double org_t_st; //原種波形のうち採用開始時刻
            double org_t_ed; //原種波形のうち採用完了時刻
            int tsp_nd;      //目標スペクトル定義点数
            
            double work1,work2;
            int kc;
            int i,j,k,kkk;
            double x1 = 0.0, y1 = 0.0, x2 = 0.0, y2 = 0.0, xx = 0.0;

            System.IO.StreamReader sr;
            System.IO.StreamWriter sw;
            string[] sbuf;
            char delim = ',';
            string dat;
            string fname1 = "";   //原種波形入力ファイル（csv）
            string fname2 = "";   //目標スペクトル入力ファイル（csv）
            string fname3 = "";   //応答スペクトル出力ファイル（csv）
            string fname4 = "";   //模擬波形時刻歴出力ファイル（csv）
            string fnameFIG = ""; //画像保存ファイル名（png）
            DRAWPARA pltdata;     //作図用パラメータ

            //*********************************************
            //データ読み込み
            //*********************************************
            //+++++++++++++++++++++++++++++++++++++++++++++
            //原種波形読み込み
            //+++++++++++++++++++++++++++++++++++++++++++++
            openFileDialog1.Title = "原種波形読み込み";
            openFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory();
            if (openFileDialog1.ShowDialog() == DialogResult.OK) fname1 = openFileDialog1.FileName;
            sr = new System.IO.StreamReader(fname1, System.Text.Encoding.Default);
            dat = sr.ReadLine();
            dat = sr.ReadLine(); sbuf = dat.Split(delim); dt = double.Parse(sbuf[1]);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); ndata = int.Parse(sbuf[1]);
            double[] org_wave = new double[ndata];//原種波形時刻歴加速度
            for (i = 1; i <= ndata; i++)
            {
                dat = sr.ReadLine(); sbuf = dat.Split(delim); org_wave[i - 1] = double.Parse(sbuf[0]);
            }
            sr.Close();
            label1.Text = fname1;
            //+++++++++++++++++++++++++++++++++++++++++++++
            //原種波形時刻歴画面描画
            //+++++++++++++++++++++++++++++++++++++++++++++
            double[] ddy = new double[ndata]; //加速度時刻歴
            double[] dy  = new double[ndata]; //速度時刻歴
            double[] y   = new double[ndata]; //変位時刻歴
            string[] scom=new string[3];
            int ncase;
            int[] nda=new int [3];
            double[,] datax=new double[3,ndata];
            double[,] datay=new double[3,ndata];

            ncase = 0;
            scom[0] = "原種波形";
            nda[0] = ndata - 1;
            ddymax = 0.0;
            for(i=0;i<=nda[0];i++){
                datax[0, i] = dt * (double)i;
                datay[0, i] = org_wave[i];
                if(ddymax<Math.Abs(org_wave[i]))ddymax = Math.Abs(org_wave[i]);
            }
            org_t_st = 0.0;
            org_t_ed = 60.0;
            DEC_PLTDATA1(out pltdata, ddymax, org_t_st, org_t_ed);
            PLTCTL(pltdata, ncase, nda, datax, datay, scom);

            //+++++++++++++++++++++++++++++++++++++++++++++
            //原種波形時刻範囲選定
            //+++++++++++++++++++++++++++++++++++++++++++++
            //原種波形中で採用する時刻の開始と終了を設定する
            int k_start, k_end;
            string setime;

            setime = Interaction.InputBox("開始時刻，終了時刻", "原種波形採用時刻","",-1,-1);
            sbuf = setime.Split(delim);
            org_t_st = double.Parse(sbuf[0]);
            org_t_ed = double.Parse(sbuf[1]);
            k_start = (int)(org_t_st / dt) + 1;
            k_end = (int)(org_t_ed / dt) + 1;
            if (ndata < k_end) k_end = ndata;
            ndata = k_end - k_start + 1;
            for (i = 1; i <= ndata; i++)
            {
                org_wave[i - 1] = org_wave[k_start + i - 2];
            }
            for (i = 0; i <= ndata-1; i++)
            {
                datax[0, i] = dt * (double)i;
                datay[0, i] = org_wave[i];
                if (ddymax < Math.Abs(org_wave[i])) ddymax = Math.Abs(org_wave[i]);
            }
            DEC_PLTDATA1(out pltdata, ddymax, 0.0, dt * (double)ndata);
            PLTCTL(pltdata, ncase, nda, datax, datay, scom);

            //+++++++++++++++++++++++++++++++++++++++++++++
            //振動の開始と終端の指定
            //+++++++++++++++++++++++++++++++++++++++++++++
            //振動開始時刻前・振動終端時刻以降は模擬波形加速度を0セットする
            double wsta, wend;
            setime = Interaction.InputBox("振動開始時刻，振動終了時刻", "振動開始と終端指定", "", -1, -1);
            sbuf = setime.Split(delim);
            wsta = double.Parse(sbuf[0]);
            wend = double.Parse(sbuf[1]);
      
            //+++++++++++++++++++++++++++++++++++++++++++++
            //目標スペクトル読み込み
            //+++++++++++++++++++++++++++++++++++++++++++++
            openFileDialog1.Title = "目標スペクトル読み込み";
            openFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory();
            if (openFileDialog1.ShowDialog() == DialogResult.OK) fname2 = openFileDialog1.FileName;
            sr = new System.IO.StreamReader(fname2, System.Text.Encoding.Default);
            dat = sr.ReadLine();
            dat = sr.ReadLine(); sbuf = dat.Split(delim); tsp_nd = int.Parse(sbuf[1]);
            dat = sr.ReadLine();
            double[] tsp_sec = new double[tsp_nd];//目標スペクトル定義周期
            double[] tsp_acc = new double[tsp_nd];//目標スペクトル定義加速度
            for (i = 1; i <= tsp_nd; i++)
            {
                dat = sr.ReadLine(); sbuf = dat.Split(delim);
                tsp_sec[i - 1] = double.Parse(sbuf[0]);
                tsp_acc[i - 1] = double.Parse(sbuf[1]);
            }
            sr.Close();
            label2.Text = fname2;
            //*********************************************
            //解析準備
            //*********************************************
            //データ数を2の累乗にセット
            nn = 2;
            do{
                nn = nn * 2;
            }while(nn < ndata * 3); //必要に応じて後続の0を追加
            //時刻歴波形初期化
            for(i=1;i<=ndata;i++){
                ddy[i - 1] = org_wave[i - 1];
            }
            //応答スペクトルを計算する固有周期(sec)の算定（周期を小さい順に格納）
            nfold = (int)(nn / 2) + 1; //加速度応答スペクトル計算周期個数

            double[] xr     =new double[nn];    //フーリエ変換・逆変換の実数部
            double[] xi     =new double[nn];    //フーリエ変換・逆変換の係数虚数部
            double[] TK     =new double[nfold]; //応答加速度を計算する周期（フーリエ変換する周波数に依存）
            double[] TSGAL  =new double[nfold]; //応答加速度を計算する周期での目標スペクトル加速度値
            double[] resacc =new double[nfold]; //応答加速度
            double[] resvel =new double[nfold]; //応答速度
            double[] resdis =new double[nfold]; //応答変位
            double[] org_res=new double[nfold]; //原種波形加速度応答格納

            for(i=nfold;i>=2;i--){
                TK[nfold - i] = dt * (double)nn / (double)(i - 1);
            }
            TK[nfold - 1] = 2.0 * dt * (double)nn;
            //目標スペクトルの加速度値セット
            for(i=0;i<=nfold-1;i++){
                if(TK[i]<tsp_sec[0]){
                    x1 = Math.Log10(tsp_sec[0]);
                    x2 = Math.Log10(tsp_sec[1]);
                    y1 = Math.Log10(tsp_acc[0]);
                    y2 = Math.Log10(tsp_acc[1]);
                }
                if(tsp_sec[0]<=TK[i]&&TK[i]<tsp_sec[tsp_nd - 1]){
                    for(j=0;j<=tsp_nd-2;j++){
                        if(tsp_sec[j]<=TK[i]&&TK[i]<tsp_sec[j + 1]){
                            x1 = Math.Log10(tsp_sec[j]);
                            x2 = Math.Log10(tsp_sec[j + 1]);
                            y1 = Math.Log10(tsp_acc[j]);
                            y2 = Math.Log10(tsp_acc[j + 1]);
                        }
                    }
                }
                if(tsp_sec[tsp_nd - 1]<=TK[i]){
                    x1 = Math.Log10(tsp_sec[tsp_nd - 2]);
                    x2 = Math.Log10(tsp_sec[tsp_nd - 1]);
                    y1 = Math.Log10(tsp_acc[tsp_nd - 2]);
                    y2 = Math.Log10(tsp_acc[tsp_nd - 1]);
                }
                xx = Math.Log10(TK[i]);
                TSGAL[i] = Math.Pow(10.0,CAL_YY(x1, y1, x2, y2, xx));
            }

            //*********************************************
            //収束計算
            //*********************************************
            for(kkk=1;kkk<=100;kkk++){
                //+++++++++++++++++++++++++++++++++++++++++++++
                //時刻歴波形の応答スペクトル計算
                //+++++++++++++++++++++++++++++++++++++++++++++
                IACC(ndata, dt, ddy, dy, y, out ddymax, out dymax, out ymax); //加速度時刻歴数値積分
                CRAC(ndata, dt, ddy, dy, y, ddymax, dymax, ymax); //加速度時刻歴基線補正
                ERES(ndata, dt, h, ddy, nfold, TK, resacc, resvel, resdis); //応答スペクトル計算
                if(kkk==1){
                    for(i=1;i<=nfold;i++){
                        org_res[i - 1] = resacc[i - 1];
                    }
                }
                //+++++++++++++++++++++++++++++++++++++++++++++
                //収束判定
                //+++++++++++++++++++++++++++++++++++++++++++++
                kc = 0;
                work1 = 0.0 ; work2 = 0.0;
                for(i=1;i<=nfold;i++){
                    rk = TSGAL[nfold - i] / resacc[nfold - i];
                    if(tsp_sec[0]<=TK[nfold - i]&&TK[nfold - i]<=tsp_sec[tsp_nd - 1]){
                        kc = kc + 1;
                        work1 = work1 + (1.0 - rk)*(1.0 - rk);
                    }
                    work2 = work2 + (1.0 - rk)*(1.0 - rk);
                }
                eps1 = Math.Sqrt(work1 / (double)kc);//目標スペクトル定義周期範囲誤差
                eps2 = Math.Sqrt(work2 / (double)nfold); //全周期範囲誤差
                if(eps1<=0.05&&eps2<= 0.05)break;
                //+++++++++++++++++++++++++++++++++++++++++++++
                //フーリエ変換
                //+++++++++++++++++++++++++++++++++++++++++++++
                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] = ddy[i - 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;
                }
                //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
                //フーリエ係数補正
                //    応答スペクトル（加速度）は周期の小さい順に格納されているが，
                //    フーリエ変換は振動数が小さい順に格納されている．
                //    このためフーリエ係数配列の添字を小→大の順とすることに対応し，
                //    応答加速度配列の添字は大→小の順として対応をとる．
                //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
                for(i=1;i<=nfold;i++){
                    rk = TSGAL[nfold - i] / resacc[nfold - i];
                    xr[i - 1] = rk * xr[i - 1];
                    xi[i - 1] = rk * xi[i - 1];
                }
                //フーリエ係数の対称性の保証
                for(i=2;i<=nfold-1;i++){
                    xr[nn - i + 1] = xr[i - 1];
                    xi[nn - i + 1] = -xi[i - 1];
                }
                xr[0] = 0.0;
                xi[0] = 0.0;
                xi[nfold - 1] = 0.0;
                //+++++++++++++++++++++++++++++++++++++++++++++
                //フーリエ逆変換
                //+++++++++++++++++++++++++++++++++++++++++++++
                for(i=1;i<=nn;i++){ //フーリエ逆変換用に虚数部の符号を反転
                    xr[i - 1] = xr[i - 1];
                    xi[i - 1] = -xi[i - 1];
                }
                FFT(nn, xr, xi); //フーリエ逆変換
                //+++++++++++++++++++++++++++++++++++++++++++++
                //時刻歴波形セット（時刻歴の最大加速度は最小周期応答加速度に一致させる）
                //+++++++++++++++++++++++++++++++++++++++++++++
                for(i=1;i<=ndata;i++){
                    ddy[i - 1] = xr[i - 1];
                    if(TSGAL[0]<Math.Abs(ddy[i - 1])){
                        if(0.0<ddy[i - 1])ddy[i - 1] = Math.Ceiling(TSGAL[0]);
                        if(ddy[i - 1]<0.0)ddy[i - 1] = -Math.Ceiling(TSGAL[0]);
                    }
                }
                //+++++++++++++++++++++++++++++++++++++++++++++
                //振動の開始と終端のセット
                //+++++++++++++++++++++++++++++++++++++++++++++
                for (i = 1; i <= ndata; i++){
                    if(dt*(double)(i-1)<wsta)ddy[i - 1] = 0.0;
                    if(wend<dt*(double)(i-1))ddy[i - 1] = 0.0;
                }
                
            }
            //*********************************************
            //収束計算完了
            //*********************************************
            toolStripLabel1.Text = "kkk=" + kkk.ToString("0");
            toolStripLabel2.Text = "eps1=" + eps1.ToString("0.000");
            toolStripLabel3.Text = "eps2=" + eps2.ToString("0.000");
            //*********************************************
            //結果ファイル出力
            //*********************************************
            fname3 = System.IO.Path.GetDirectoryName(label2.Text) + "\\mogi-spec.csv";
            fname4 = System.IO.Path.GetDirectoryName(label2.Text) + "\\mogi-wave.csv";
            //+++++++++++++++++++++++++++++++++++++++++++++
            //模擬波応答スペクトル出力
            //+++++++++++++++++++++++++++++++++++++++++++++
            sw = new System.IO.StreamWriter(fname3, false, System.Text.Encoding.Default);
            sw.WriteLine("模擬波応答スペクトル");
            sw.WriteLine("ndata," + nfold.ToString("0"));
            sw.WriteLine("Period(sec),acc(gal)");
            for(i=1;i<=nfold;i++){
                dat = TK[i - 1].ToString("0.000000");
                dat = dat + "," + resacc[i - 1].ToString("0.0");
                sw.WriteLine(dat);
            }
            sw.Close();
            //+++++++++++++++++++++++++++++++++++++++++++++
            //模擬波時刻歴出力
            //+++++++++++++++++++++++++++++++++++++++++++++
            sw = new System.IO.StreamWriter(fname4, false, System.Text.Encoding.Default);
            sw.WriteLine("模擬波時刻歴");
            sw.WriteLine("dt," + dt.ToString("0.00"));
            sw.WriteLine("ndata," + ndata.ToString("0"));
            for(i=1;i<=ndata;i++){
                dat = ddy[i - 1].ToString("0.0");
                sw.WriteLine(dat);
            }
            sw.Close();
            //+++++++++++++++++++++++++++++++++++++++++++++
            //原種波形時刻歴画像出力
            //+++++++++++++++++++++++++++++++++++++++++++++
            ncase = 0; 
            k = 0;
            scom[k] = "原種波形（" + org_t_st.ToString("0") + "〜" + org_t_ed.ToString("0") + " sec）";
            nda[k] = ndata - 1;
            ddymax = 0.0;
            for(i=0;i<=nda[k];i++){
                datax[k,i] = dt * (double)i;
                datay[k,i] = org_wave[i];
                if(ddymax<Math.Abs(org_wave[i]))ddymax = Math.Abs(org_wave[i]);
            }
            DEC_PLTDATA1(out pltdata, ddymax, org_t_st, org_t_ed);
            PLTCTL(pltdata, ncase, nda, datax, datay, scom);
            fnameFIG = System.IO.Path.GetDirectoryName(label2.Text) + "\\fig-0-wave.png";
            pictureBox2.Image.Save(fnameFIG, System.Drawing.Imaging.ImageFormat.Png);
            //+++++++++++++++++++++++++++++++++++++++++++++
            //模擬波時刻歴画像出力
            //+++++++++++++++++++++++++++++++++++++++++++++
            ncase = 0;
            k=0;
            scom[k] = "模擬波形";
            nda[k] = ndata - 1;
            ddymax = 0.0;
            for(i=0;i<=nda[k];i++){
                datax[k,i] = dt * (double)i;
                datay[k,i] = ddy[i];
                if(ddymax<Math.Abs(ddy[i]))ddymax = Math.Abs(ddy[i]);
            }
            DEC_PLTDATA1(out pltdata, ddymax, org_t_st, org_t_ed);
            PLTCTL(pltdata, ncase, nda, datax, datay, scom);
            fnameFIG = System.IO.Path.GetDirectoryName(label2.Text) + "\\fig-1-wave.png";
            pictureBox2.Image.Save(fnameFIG, System.Drawing.Imaging.ImageFormat.Png);
            //+++++++++++++++++++++++++++++++++++++++++++++
            //加速度応答スペクトル
            //+++++++++++++++++++++++++++++++++++++++++++++
            ncase = 2;
            k = 0;
            scom[k] = "原種波形";
            nda[k] = nfold - 1;
            for(i=0;i<=nda[k];i++){
                datax[k, i] = TK[i];
                datay[k, i] = org_res[i];
            }
            k = 1;
            scom[k] = "目標SP";
            nda[k] = tsp_nd - 1;
            for(i=0;i<=nda[k];i++){
                datax[k, i] = tsp_sec[i];
                datay[k, i] = tsp_acc[i];
            }
            k = 2;
            scom[k] = "模擬波形";
            nda[k] = nfold - 1;
            for(i=0;i<=nda[k];i++){
                datax[k, i] = TK[i];
                datay[k, i] = resacc[i];
            }
            DEC_PLTDATA2(out pltdata);
            PLTCTL(pltdata, ncase, nda, datax, datay, scom);
            fnameFIG = System.IO.Path.GetDirectoryName(label2.Text) + "\\fig-2-sp.png";
            pictureBox2.Image.Save(fnameFIG, System.Drawing.Imaging.ImageFormat.Png);
        }
        //---------------------------------------------------------------------------------
        private double CAL_YY(double x1, double y1, double x2, double y2, double xx)
        {
            //２点を通る直線による内挿
            double aa, bb;
            aa = (y2 - y1) / (x2 - x1);
            bb = (x2 * y1 - x1 * y2) / (x2 - x1);
            return aa * xx + bb;
        }
        //---------------------------------------------------------------------------------
        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 ERES(int nn, double dt, double h, double[] ddy, int nt,
                          double[] tk, double[] resacc, double[] resvel, double[] resdis)
        {
            //応答値算出
            //nn      :加速度時刻歴総数
            //dt      :時間間隔（入力値）
            //h       :減衰定数（入力値）
            //ddy[]   :加速度時刻歴（入力値）
            //nt      :計算周期総数[データ格納0〜nt-1]（入力値）
            //tk[]    :計算周期（入力値）
            //resacc[]:加速度応答スペクトル（計算値）
            //resvel[]:速度応答スペクトル（計算値）
            //resdis[]:変位応答スペクトル（計算値）

            int i, m;
            double accmax, velmax, dismax;//応答最大値(gal,kine,cm)
            double w, w2, hw, wd, wdt, e, cwdt, swdt;
            double ss, cc, s1, c1, s2, c2, s3, c3;
            double A11, A12, A21, A22, B11, B12, B21, B22;
            double dxf, xf, ddym, ddyf, x, dx, ddx;

            for (i = 0; i <= nt - 1; i++)
            {
                w = 2.0 * pi / tk[i];
                w2 = w * w;
                hw = h * w;
                wd = w * Math.Sqrt(1.0 - h * h);
                wdt = wd * dt;
                e = Math.Exp(-hw * dt);
                cwdt = Math.Cos(wdt);
                swdt = Math.Sin(wdt);
                A11 = e * (cwdt + hw * swdt / wd);
                A12 = e * swdt / wd;
                A21 = -e * w2 * swdt / wd;
                A22 = e * (cwdt - hw * swdt / wd);
                ss = -hw * swdt - wd * cwdt;
                cc = -hw * cwdt + wd * swdt;
                s1 = (e * ss + wd) / w2;
                c1 = (e * cc + hw) / w2;
                s2 = (e * dt * ss + hw * s1 + wd * c1) / w2;
                c2 = (e * dt * cc + hw * c1 - wd * s1) / w2;
                s3 = dt * s1 - s2;
                c3 = dt * c1 - c2;
                B11 = -s2 / wdt;
                B12 = -s3 / wdt;
                B21 = (hw * s2 - wd * c2) / wdt;
                B22 = (hw * s3 - wd * c3) / wdt;
                accmax = 2.0 * hw * Math.Abs(ddy[0]) * dt;
                velmax = Math.Abs(ddy[0]) * dt;
                dismax = 0.0;
                dxf = -ddy[0] * dt;
                xf = 0.0;
                for (m = 1; m <= nn - 1; m++)
                {
                    ddym = ddy[m];
                    ddyf = ddy[m - 1];
                    x = A12 * dxf + A11 * xf + B12 * ddym + B11 * ddyf;
                    dx = A22 * dxf + A21 * xf + B22 * ddym + B21 * ddyf;
                    ddx = -2.0 * hw * dx - w2 * x;
                    if (accmax <= Math.Abs(ddx)) accmax = Math.Abs(ddx);
                    if (velmax <= Math.Abs(dx)) velmax = Math.Abs(dx);
                    if (dismax <= Math.Abs(x)) dismax = Math.Abs(x);
                    dxf = dx;
                    xf = x;
                }
                resacc[i] = accmax;
                resvel[i] = velmax;
                resdis[i] = dismax;
            }
        }
        //---------------------------------------------------------------------------
        private void IACC(int nn, double dt, double[] ddy, double[] dy, double[] y,
                          out double ddymax, out double dymax, out double ymax)
        {
            //加速度時刻歴数値積分
            //nn:加速度時刻歴データ総数
            //dt    :時間間隔（入力値）
            //ddy[] :加速度時刻歴（入力値）
            //dy[]  :速度時刻歴（計算出力）
            //y[]   :変位時刻歴（計算出力値）
            //ddymax:加速度最大値（計算出力値）
            //dymax :速度最大値（計算出力値）
            //ymax  :変位最大値（計算出力値）

            int i;

            ddymax = 0.0;
            dymax = 0.0;
            ymax = 0.0;
            dy[0] = 0.0;
            y[0] = 0.0;
            for (i = 1; i <= nn - 1; i++)
            {
                dy[i] = dy[i - 1] + (ddy[i - 1] + ddy[i]) * dt / 2.0;
                y[i] = y[i - 1] + dy[i - 1] * dt + (ddy[i - 1] / 3.0 + ddy[i] / 6.0) * dt * dt;
                if (ddymax < Math.Abs(ddy[i])) ddymax = Math.Abs(ddy[i]);
                if (dymax < Math.Abs(dy[i])) dymax = Math.Abs(dy[i]);
                if (ymax < Math.Abs(y[i])) ymax = Math.Abs(y[i]);
            }
        }
        //---------------------------------------------------------------------------
        private void CRAC(int nn, double dt, double[] ddy, double[] dy, double[] y,
                          double ddymax, double dymax, double ymax)
        {
            //加速度時刻歴基線補正
            //nn:加速度時刻歴データ総数
            //dt    :時間間隔（入力値）
            //ddy[] :加速度時刻歴（入力値・書換出力値）
            //dy[]  :速度時刻歴（入力値）
            //y[]   :変位時刻歴（入力値）
            //ddymax:加速度最大値（入力値）
            //dymax :速度最大値（入力値）
            //ymax  :変位最大値（入力値）

            int i;
            double tt, t, sum, a1, a0, acmax, coef;

            tt = (double)(nn - 1) * dt;
            t = 0.0;
            for (i = 0; i <= nn - 1; i++)
            {
                y[i] = y[i] * (3.0 * tt - 2.0 * t) * t * t;
                t = t + dt;
            }
            sum = (y[0] + y[nn - 1]) / 2.0;
            for (i = 1; i <= nn - 2; i++)
            {
                sum = sum + y[i];
            }
            sum = sum * dt;
            a1 = 28.0 / 13.0 / tt / tt * (2.0 * dy[nn - 1] - 15.0 / Math.Pow(tt, 5.0) * sum);
            a0 = dy[nn - 1] / tt - a1 / 2.0 * tt;
            t = 0.0;
            acmax = 0.0;
            for (i = 0; i <= nn - 1; i++)
            {
                ddy[i] = ddy[i] - a0 - a1 * t;
                if (acmax < Math.Abs(ddy[i])) acmax = Math.Abs(ddy[i]);
                t = t + dt;
            }
            coef = ddymax / acmax;
            for (i = 0; i <= nn - 1; i++)
            {
                ddy[i] = ddy[i] * coef;
            }
        }
        //---------------------------------------------------------------------------
        private void PLTCTL(DRAWPARA pltdata, int ncase,int[] nda,
                            double[,] datax, double[,] datay, string[] scom)
        {
            //nda[]：データ数（入力はデータ数，プログラム上数値はデータ数−１）
            //datax[,]：プロットデータｘ座標
            //datay[,]：プロットデータｙ座標
            int i,kk;
            PictureBox thePicBox = pictureBox1;
            Bitmap bmp;
            Graphics g;

            //データ対数値処理
            for (kk = 0; kk <= ncase;kk++){
                if (pltdata.slogx == "L")
                {
                    for (i = 0; i <= nda[kk]; i++)
                    {
                        if (datax[kk, i] < Math.Pow(10.0, pltdata.xmin)) datax[kk, i] = Math.Pow(10.0, pltdata.xmin);
                        datax[kk, i] = Math.Log10(datax[kk, i]);
                    }
                }
                if (pltdata.slogy == "L")
                {
                    for (i = 0; i <= nda[kk]; i++)
                    {
                        if (datay[kk,i] < Math.Pow(10.0, pltdata.ymin))datay[kk,i] = Math.Pow(10.0, pltdata.ymin); 
                        datay[kk,i] = Math.Log10(datay[kk,i]);
                    }
                }
            }

            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, ncase, nda, datax, datay, scom);
                g.Dispose();
            }
            pictureBox1.Left = 0;
            pictureBox1.Top = toolStrip1.Height;
            this.ClientSize = new Size(pictureBox1.Width, pictureBox1.Height + toolStrip1.Height);

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

            int ds=40;
            int twmax=0;

            //凡例描画余白確保
            if(0<ncase){
                twmax = 0;
                for(i=0;i<=ncase;i++){
                    TextSize1 = g.MeasureString(scom[i], f);
                    if(twmax<(int)TextSize1.Width)twmax = (int)TextSize1.Width;
                }
                twmax = twmax + (ds + 20) * pltdata.kpt;
            }
            kxi = pltdata.kpt * pltdata.kxi0;
            kxf = pltdata.kpt * pltdata.kxf0 - twmax;
            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;
            if (pltdata.slogx == "N")DRN_XJIKU(g, f, LPen, pltdata, kxi, kxf, kyi, kyf);
            if (pltdata.slogy == "N")DRN_YJIKU(g, f, LPen, pltdata, kxi, kxf, kyi, kyf);
            if (pltdata.slogx == "L")DRL_XJIKU(g, f, LPen, pltdata, kxi, kxf, kyi, kyf);
            if (pltdata.slogy == "L")DRL_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);
            //データを線で連結
            LPen.DashStyle = System.Drawing.Drawing2D.DashStyle.Solid;
            LPen.Color = Color.Blue; 
            LPen.Width = 1;
            for (kk = 0; kk <= ncase; kk++)
            {
                if (ncase != 0)
                {
                    switch (kk)
                    {
                        case 0: LPen.Color = Color.Lime; LPen.Width = 1; break;
                        case 1: LPen.Color = Color.Red; LPen.Width = 5; break;
                        case 2: LPen.Color = Color.Blue; LPen.Width = 2; break;
                    }
                }
                xx = datax[kk,0];
                yy = datay[kk,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 <= nda[kk]; i++)
                {
                    xx = datax[kk,i];
                    yy = datay[kk,i];
                    if (pltdata.xmin <= xx && xx <= pltdata.xmax)
                    {

                        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(LPen, kxx1, kyy1, kxx2, kyy2);
                    }
                    kxx1 = kxx2;
                    kyy1 = kyy2;
                }
            }
            //コメント・凡例描画
            if (ncase == 0){
                //右上コメント
                kk = 0;
                str = scom[kk];
                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);
            }else{
                //右側凡例
                for(i=0;i<=ncase;i++){
                    str =scom[i];
                    TextSize1 = g.MeasureString(str, f);
                    g.DrawString(str, f, Brushes.Black, kxf + pltdata.kpt * (ds + 20), kyi + TextSize1.Height * i);
                    switch(i){
                        case 0 : LPen.Color = Color.Lime ; LPen.Width = 1; break;
                        case 1 : LPen.Color = Color.Red  ; LPen.Width = 5; break;
                        case 2 : LPen.Color = Color.Blue ; LPen.Width = 2; break;
                    }
                    kxx1 = kxf + pltdata.kpt * 10;
                    kxx2 = kxx1 + ds * pltdata.kpt;
                    kyy1 = kyi + (int)(TextSize1.Height * i) + (int)(TextSize1.Height / 2);
                    kyy2 = kyy1;
                    g.DrawLine(LPen, kxx1, kyy1, kxx2, kyy2);
                }
                str = "PGA（gal）";
                TextSize1 = g.MeasureString(str, f);
                kyy1 = kyi + (int)(TextSize1.Height * 5);
                g.DrawString(str, f, Brushes.Black, kxf + pltdata.kpt * 10, kyy1);
                for(i=0;i<=ncase;i++){
                    str = Math.Pow(10,datay[i, 0]).ToString("0") + "：" + scom[i];
                    TextSize1 = g.MeasureString(str, f);
                    g.DrawString(str, f, Brushes.Black, kxf + pltdata.kpt * 10, kyy1 + TextSize1.Height * (i + 1));
                }
            }
            f.Dispose();
            LPen.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 DRL_XJIKU(Graphics g, System.Drawing.Font f, System.Drawing.Pen LPen,
                               DRAWPARA pltdata, int kxi, int kxf, int kyi, int kyf)
        {
            int i, j, kxx;
            double xx;
            string str;
            System.Drawing.SizeF TextSize1;
            //対数ｘ軸描画
            for (i = (int)pltdata.xmin; i <= (int)pltdata.xmax - 1; i++)
            {
                for (j = 1; j <= 9; j++)
                {
                    xx = Math.Log10((double)(j) * Math.Pow(10.0, (double)i));
                    kxx = kxi + (int)((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin));
                    g.DrawLine(LPen, kxx, kyi, kxx, kyf);
                }
            }
            for (i = (int)pltdata.xmin; i <= (int)pltdata.xmax; i++)
            {
                xx = (double)i;
                kxx = kxi + (int)((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin));
                str = Math.Pow(10.0, (double)i).ToString();
                TextSize1 = g.MeasureString(str, f);
                g.DrawString(str, f, Brushes.Black, kxx - TextSize1.Width / 2, kyf + pltdata.kpt * 3);
            }
        }
        //---------------------------------------------------------------------------------
        private void DRL_YJIKU(Graphics g, System.Drawing.Font f, System.Drawing.Pen LPen,
                               DRAWPARA pltdata, int kxi, int kxf, int kyi, int kyf)
        {
            int i, j;
            double yy;
            string str;
            int kyy;
            System.Drawing.SizeF TextSize1;
            //対数ｙ軸描画
            for (i = (int)pltdata.ymin; i <= (int)pltdata.ymax - 1; i++)
            {
                for (j = 1; j <= 9; j++)
                {
                    yy = Math.Log10((double)j * Math.Pow(10.0, (double)i));
                    kyy = kyf - (int)((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin));
                    g.DrawLine(LPen, kxi, kyy, kxf, kyy);
                }
            }
            for (i = (int)pltdata.ymin; i <= (int)pltdata.ymax; i++)
            {
                yy = (double)i;
                kyy = kyf - (int)((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin));
                str = (Math.Pow(10.0, (double)i)).ToString();
                TextSize1 = g.MeasureString(str, f);
                g.DrawString(str, f, Brushes.Black, kxi - TextSize1.Width - pltdata.kpt * 2, kyy - TextSize1.Height / 2);
            }
        }
        //---------------------------------------------------------------------------------
        private void DEC_PLTDATA1(out DRAWPARA pltdata, double ddymax,
                         double org_t_st, double org_t_ed)
        {
            //波形時刻歴用作図用データ
            pltdata.kpt = 1;
            pltdata.HSIZE = 600;
            pltdata.DHL = 70;
            pltdata.DHR = 10;
            pltdata.VSIZE = 150;
            pltdata.DVL = 40;
            pltdata.DVU = 20;
            pltdata.sxjiku = "時刻（sec）";
            pltdata.xmin = 0.0;
            pltdata.xmax = org_t_ed - org_t_st;
            pltdata.dx = 5.0;
            pltdata.slogx = "N";
            pltdata.syjiku = "加速度（gal）";
            pltdata.ymin = -Math.Ceiling(ddymax);
            pltdata.ymax = Math.Ceiling(ddymax);
            pltdata.dy = Math.Ceiling(ddymax);
            pltdata.slogy = "N";
            //プロット領域定義
            pltdata.kxi0 = pltdata.DHL;
            pltdata.kxf0 = pltdata.HSIZE - pltdata.DHR;
            pltdata.kyi0 = pltdata.DVU;
            pltdata.kyf0 = pltdata.VSIZE - pltdata.DVL;
        }
        //---------------------------------------------------------------------------------
        private void DEC_PLTDATA2(out DRAWPARA pltdata)
        {
            //加速度応答スペクトル作図用データ
            pltdata.kpt = 1;
            pltdata.HSIZE = 600;
            pltdata.DHL = 70;
            pltdata.DHR = 10;
            pltdata.VSIZE = 300;
            pltdata.DVL = 40;
            pltdata.DVU = 20;
            pltdata.sxjiku = "周期（sec）";
            pltdata.xmin = -2.0;
            pltdata.xmax = 1.0;
            pltdata.dx = 1.0;
            pltdata.slogx = "L";
            pltdata.syjiku = "加速度応答スペクトル（gal）";
            pltdata.ymin = 1.0;
            pltdata.ymax = 4.0;
            pltdata.dy = 1.0;
            pltdata.slogy = "L";
            //プロット領域定義
            pltdata.kxi0 = pltdata.DHL;
            pltdata.kxf0 = pltdata.HSIZE - pltdata.DHR;
            pltdata.kyi0 = pltdata.DVU;
            pltdata.kyf0 = pltdata.VSIZE - pltdata.DVL;
        }
        //---------------------------------------------------------------------------------

    }
}