Option Explicit On
Option Strict On
Public Class Form1

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

    Private Structure DRAWVALUE
        Dim ncase As Integer    'プロットする線の数から１を減じた数値
        Dim nda() As Integer    '１曲線を描くデータ組数から１を減じた数値
        Dim datax(,) As Double  '曲線を描くデータの横軸値
        Dim datay(,) As Double  '曲線を描くデータの縦軸値
        Dim scom() As String    '曲線の名称（凡例に使用）
    End Structure


    Private Sub Form1_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load
        ToolStripLabel1.Text = "繰返回数"
        ToolStripLabel2.Text = "誤差（目標SP定義範囲）"
        ToolStripLabel3.Text = "誤差（全体）"
    End Sub


    Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click
        ToolStripLabel1.Text = "繰返回数"
        ToolStripLabel2.Text = "誤差（目標SP定義範囲）"
        ToolStripLabel3.Text = "誤差（全体）"
        Call main()
        MessageBox.Show("計算完了", "完了通知")
    End Sub


    Private Function CAL_YY(ByVal x1 As Double, ByVal y1 As Double, _
                            ByVal x2 As Double, ByVal y2 As Double, _
                            ByRef xx As Double) As Double
        '２点を通る直線による内挿
        Dim aa As Double
        Dim bb As Double
        aa = (y2 - y1) / (x2 - x1)
        bb = (x2 * y1 - x1 * y2) / (x2 - x1)
        CAL_YY = aa * xx + bb
    End Function


    Private Sub DEC_PLTDATA1(ByRef pltdata As DRAWPARA, ByVal ddymax As Double, _
                             ByVal org_t_st As Double, ByVal org_t_ed As Double)
        '波形時刻歴用作図用データ
        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"
    End Sub


    Private Sub DEC_PLTDATA2(ByRef pltdata As DRAWPARA)
        '加速度応答スペクトル作図用データ
        pltdata.kpt = 1
        pltdata.HSIZE = 600
        pltdata.DHL = 70
        pltdata.DHR = 10
        pltdata.VSIZE = 300
        pltdata.DVL = 40
        pltdata.DVU = 10
        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"
    End Sub


    Private Sub main()

        Dim dt As Double '時刻歴波形サンプリングピッチ
        Dim ndata As Integer '時刻歴波形データ数
        Dim org_wave() As Double '原種波形時刻歴加速度
        Dim org_t_st As Double '原種波形のうち採用開始時刻
        Dim org_t_ed As Double '原種波形のうち採用完了時刻
        Dim tsp_nd As Integer '目標スペクトル定義点数
        Dim tsp_sec() As Double '目標スペクトル定義周期
        Dim tsp_acc() As Double '目標スペクトル定義加速度
        Dim nn As Integer 'FFTにかけるデータ数（２の累乗）
        Dim xr() As Double 'フーリエ変換・逆変換の実数部
        Dim xi() As Double 'フーリエ変換・逆変換の係数虚数部
        Dim TK() As Double '応答加速度を計算する周期（フーリエ変換する周波数に依存）
        Dim TSGAL() As Double '応答加速度を計算する周期での目標スペクトル加速度値
        Dim nfold As Integer '折り曲げ振動数を与える番号に１を加算した数値

        Dim h As Double = 0.05 '減衰定数
        Dim ddy() As Double '加速度時刻歴
        Dim dy() As Double '速度時刻歴
        Dim y() As Double '変位時刻歴
        Dim resacc() As Double '応答加速度
        Dim resvel() As Double '応答速度
        Dim resdis() As Double '応答変位
        Dim ddymax As Double '最大加速度
        Dim dymax As Double '最大速度
        Dim ymax As Double '最大変位
        Dim org_res() As Double '原種波形加速度応答格納

        Dim eps1 As Double '目標スペクトル定義周期範囲での誤差
        Dim eps2 As Double '全誤差
        Dim rk As Double '目標スペクトルと模擬地震波応答スペクトルの加速度比
        Dim work1 As Double
        Dim work2 As Double
        Dim kc As Integer

        Dim x1 As Double
        Dim y1 As Double
        Dim x2 As Double
        Dim y2 As Double
        Dim xx As Double

        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim kkk As Integer

        Dim sr As System.IO.StreamReader
        Dim sw As System.IO.StreamWriter
        Dim sbuf() As String
        Dim delim() As Char = {","c}
        Dim dat As String
        Dim fname1 As String = ""   '原種波形入力ファイル（csv）
        Dim fname2 As String = ""   '目標スペクトル入力ファイル（csv）
        Dim fname3 As String = ""   '応答スペクトル出力ファイル（csv）
        Dim fname4 As String = ""   '模擬波形時刻歴出力ファイル（csv）
        Dim fnameFIG As String = "" '画像保存ファイル名（png）

        Dim pltdata As DRAWPARA '作図用パラメータ
        Dim pltvalue As DRAWVALUE 'プロットデータ数値

        '波形時刻歴用作図用データ初期化
        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.slogx = ""
        pltdata.syjiku = ""
        pltdata.ymin = 0.0
        pltdata.ymax = 0.0
        pltdata.dy = 0.0
        pltdata.slogy = ""

        '*********************************************
        'データ読み込み
        '*********************************************
        '原種波形読み込み
        OpenFileDialog1.Title = "原種波形読み込み"
        OpenFileDialog1.InitialDirectory = My.Computer.FileSystem.CurrentDirectory
        If OpenFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fname1 = OpenFileDialog1.FileName
        sr = New System.IO.StreamReader(fname1, System.Text.Encoding.Default)
        dat = sr.ReadLine
        dat = sr.ReadLine : sbuf = dat.Split(delim) : dt = CDbl(sbuf(1))
        dat = sr.ReadLine : sbuf = dat.Split(delim) : ndata = CInt(sbuf(1))
        ReDim org_wave(ndata - 1)
        For i = 1 To ndata
            dat = sr.ReadLine : sbuf = dat.Split(delim) : org_wave(i - 1) = CDbl(sbuf(0))
        Next i
        sr.Close()
        Label1.Text = fname1
        '原種波形時刻歴画面描画
        pltvalue.ncase = 0
        ReDim pltvalue.nda(pltvalue.ncase)
        ReDim pltvalue.datax(pltvalue.ncase, ndata - 1)
        ReDim pltvalue.datay(pltvalue.ncase, ndata - 1)
        ReDim pltvalue.scom(pltvalue.ncase)
        pltvalue.scom(0) = "原種波形"
        pltvalue.nda(0) = ndata - 1
        ddymax = 0.0
        For i = 0 To pltvalue.nda(0)
            pltvalue.datax(0, i) = dt * CDbl(i)
            pltvalue.datay(0, i) = org_wave(i)
            If ddymax < Math.Abs(org_wave(i)) Then ddymax = Math.Abs(org_wave(i))
        Next i
        org_t_st = 0.0
        org_t_ed = 60.0
        Call DEC_PLTDATA1(pltdata, ddymax, org_t_st, org_t_ed)
        Call SAKUZU(pltdata, pltvalue)
        '原種波形時刻範囲選定
        '原種波形中で採用する時刻の開始と終了を設定する
        Dim k_start As Integer
        Dim k_end As Integer
        Dim setime As String
        setime = InputBox("開始時刻，終了時刻", "原種波形採用時刻")
        sbuf = setime.Split(delim)
        org_t_st = CDbl(sbuf(0))
        org_t_ed = CDbl(sbuf(1))
        k_start = CInt(org_t_st / dt) + 1
        k_end = CInt(org_t_ed / dt) + 1
        If ndata < k_end Then k_end = ndata
        ndata = k_end - k_start + 1
        For i = 1 To ndata
            org_wave(i - 1) = org_wave(k_start + i - 2)
        Next i
        '振動開始と振動終端時刻の指定
        ddymax = 0.0
        For i = 0 To ndata - 1
            pltvalue.datax(0, i) = dt * CDbl(i)
            pltvalue.datay(0, i) = org_wave(i)
            If ddymax < Math.Abs(org_wave(i)) Then ddymax = Math.Abs(org_wave(i))
        Next i
        Call DEC_PLTDATA1(pltdata, ddymax, 0.0, dt * CDbl(ndata))
        Call SAKUZU(pltdata, pltvalue)
        '振動開始時刻前と振動終端時刻以降は強制的に模擬波形加速度を0セットする
        Dim w_sta As Double '振動開始時刻
        Dim w_end As Double '振動終端時刻
        setime = InputBox("振動開始時刻，振動終了時刻", "振動開始・終了時刻設定")
        sbuf = setime.Split(delim)
        w_sta = CDbl(sbuf(0))
        w_end = CDbl(sbuf(1))
        '目標スペクトル読み込み
        OpenFileDialog1.Title = "目標スペクトル読み込み"
        OpenFileDialog1.InitialDirectory = My.Computer.FileSystem.CurrentDirectory
        If OpenFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then 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 = CInt(sbuf(1))
        dat = sr.ReadLine
        ReDim tsp_sec(tsp_nd - 1)
        ReDim tsp_acc(tsp_nd - 1)
        For i = 1 To tsp_nd
            dat = sr.ReadLine : sbuf = dat.Split(delim)
            tsp_sec(i - 1) = CDbl(sbuf(0))
            tsp_acc(i - 1) = CDbl(sbuf(1))
        Next i
        sr.Close()
        Label2.Text = fname2
        '*********************************************
        '解析準備
        '*********************************************
        'データ数を2の累乗にセット
        nn = 2
        Do
            nn = nn * 2
        Loop While nn < ndata * 3 '必要に応じて後続の0を追加
        ReDim xr(nn - 1)
        ReDim xi(nn - 1)
        ReDim ddy(ndata - 1)
        ReDim dy(ndata - 1)
        ReDim y(ndata - 1)
        '時刻歴波形初期化
        For i = 1 To ndata
            ddy(i - 1) = org_wave(i - 1)
        Next i
        '応答スペクトルを計算する固有周期(sec)の算定（周期を小さい順に格納）
        nfold = CInt(nn / 2) + 1 '加速度応答スペクトル計算周期個数
        ReDim TK(nfold - 1)
        ReDim TSGAL(nfold - 1)
        ReDim resacc(nfold - 1)
        ReDim resvel(nfold - 1)
        ReDim resdis(nfold - 1)
        ReDim org_res(nfold - 1)
        For i = nfold To 2 Step -1
            TK(nfold - i) = dt * CDbl(nn) / CDbl(i - 1)
        Next i
        TK(nfold - 1) = 2.0 * dt * CDbl(nn)
        '目標スペクトルの加速度値セット
        For i = 0 To nfold - 1
            If TK(i) < tsp_sec(0) Then
                x1 = Math.Log10(tsp_sec(0))
                x2 = Math.Log10(tsp_sec(1))
                y1 = Math.Log10(tsp_acc(0))
                y2 = Math.Log10(tsp_acc(1))
            ElseIf tsp_sec(0) <= TK(i) And TK(i) < tsp_sec(tsp_nd - 1) Then
                For j = 0 To tsp_nd - 2
                    If tsp_sec(j) <= TK(i) And TK(i) < tsp_sec(j + 1) Then
                        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))
                    End If
                Next j
            ElseIf tsp_sec(tsp_nd - 1) <= TK(i) Then
                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))
            End If
            xx = Math.Log10(TK(i))
            TSGAL(i) = 10.0 ^ (CAL_YY(x1, y1, x2, y2, xx))
        Next i

        '*********************************************
        '収束計算
        '*********************************************
        For kkk = 1 To 100
            '+++++++++++++++++++++++++++++++++++++++++++++
            '時刻歴波形の応答スペクトル計算
            '+++++++++++++++++++++++++++++++++++++++++++++
            Call IACC(ndata, dt, ddy, dy, y, ddymax, dymax, ymax) '加速度時刻歴数値積分
            Call CRAC(ndata, dt, ddy, dy, y, ddymax, dymax, ymax) '加速度時刻歴基線補正
            Call ERES(ndata, dt, h, ddy, nfold, TK, resacc, resvel, resdis) '応答スペクトル計算
            If kkk = 1 Then
                For i = 1 To nfold
                    org_res(i - 1) = resacc(i - 1)
                Next i
            End If
            '+++++++++++++++++++++++++++++++++++++++++++++
            '収束判定
            '+++++++++++++++++++++++++++++++++++++++++++++
            kc = 0
            work1 = 0.0 : work2 = 0.0
            For i = 1 To nfold
                rk = TSGAL(nfold - i) / resacc(nfold - i)
                If tsp_sec(0) <= TK(nfold - i) And TK(nfold - i) <= tsp_sec(tsp_nd - 1) Then
                    kc = kc + 1
                    work1 = work1 + (1.0 - rk) ^ 2
                End If
                work2 = work2 + (1.0 - rk) ^ 2
            Next i
            eps1 = Math.Sqrt(work1 / CDbl(kc)) '目標スペクトル定義周期範囲誤差
            eps2 = Math.Sqrt(work2 / CDbl(nfold)) '全周期範囲誤差
            If eps1 <= 0.05 And eps2 <= 0.05 Then Exit For
            '+++++++++++++++++++++++++++++++++++++++++++++
            'フーリエ変換
            '+++++++++++++++++++++++++++++++++++++++++++++
            For i = 1 To nn
                xr(i - 1) = 0.0 : xi(i - 1) = 0.0
            Next
            For i = 1 To ndata
                xr(i - 1) = ddy(i - 1)
            Next i
            Call FFT(nn, xr, xi) 'フーリエ変換
            For i = 1 To nn '戻り値をデータ数で除す
                xr(i - 1) = xr(i - 1) / CDbl(nn)
                xi(i - 1) = xi(i - 1) / CDbl(nn)
            Next i
            '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            'フーリエ係数補正
            '    応答スペクトル（加速度）は周期の小さい順に格納されているが，
            '    フーリエ変換は振動数が小さい順に格納されている．
            '    このためフーリエ係数配列の添字を小→大の順とすることに対応し，
            '    応答加速度配列の添字は大→小の順として対応をとる．
            '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            For i = 1 To nfold
                rk = TSGAL(nfold - i) / resacc(nfold - i)
                xr(i - 1) = rk * xr(i - 1)
                xi(i - 1) = rk * xi(i - 1)
            Next i
            'フーリエ係数の対称性の保証
            For i = 2 To nfold - 1
                xr(nn - i + 1) = xr(i - 1)
                xi(nn - i + 1) = -xi(i - 1)
            Next i
            xr(0) = 0.0
            xi(0) = 0.0
            xi(nfold - 1) = 0.0
            '+++++++++++++++++++++++++++++++++++++++++++++
            'フーリエ逆変換
            '+++++++++++++++++++++++++++++++++++++++++++++
            For i = 1 To nn 'フーリエ逆変換用に虚数部の符号を反転
                xr(i - 1) = xr(i - 1)
                xi(i - 1) = -xi(i - 1)
            Next i
            Call FFT(nn, xr, xi) 'フーリエ逆変換
            '+++++++++++++++++++++++++++++++++++++++++++++
            '時刻歴波形セット（時刻歴の最大加速度は最小周期応答加速度に一致させる）
            '+++++++++++++++++++++++++++++++++++++++++++++
            For i = 1 To ndata
                ddy(i - 1) = xr(i - 1)
                If TSGAL(0) < Math.Abs(ddy(i - 1)) Then
                    If 0.0 < ddy(i - 1) Then ddy(i - 1) = Math.Ceiling(TSGAL(0))
                    If ddy(i - 1) < 0.0 Then ddy(i - 1) = -Math.Ceiling(TSGAL(0))
                End If
            Next i
            '+++++++++++++++++++++++++++++++++++++++++++++
            '振動開始前・振動終了後の加速度0セット
            '+++++++++++++++++++++++++++++++++++++++++++++
            For i = 1 To ndata
                If dt * CDbl(i - 1) < w_sta Then ddy(i - 1) = 0.0
                If w_end < dt * CDbl(i - 1) Then ddy(i - 1) = 0.0
            Next i
        Next kkk
        '*********************************************
        '収束計算完了
        '*********************************************
        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 To nfold
            dat = TK(i - 1).ToString("0.000000")
            dat = dat & "," & resacc(i - 1).ToString("0.0")
            sw.WriteLine(dat)
        Next i
        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 To ndata
            dat = ddy(i - 1).ToString("0.0")
            sw.WriteLine(dat)
        Next i
        sw.Close()
        '+++++++++++++++++++++++++++++++++++++++++++++
        '原種波形時刻歴画像出力
        '+++++++++++++++++++++++++++++++++++++++++++++
        pltvalue.ncase = 0
        ReDim pltvalue.nda(pltvalue.ncase)
        ReDim pltvalue.datax(pltvalue.ncase, ndata - 1)
        ReDim pltvalue.datay(pltvalue.ncase, ndata - 1)
        ReDim pltvalue.scom(pltvalue.ncase)
        pltvalue.scom(0) = "原種波形（" & org_t_st.ToString("0") & "〜" & org_t_ed.ToString("0") & " sec）"
        pltvalue.nda(0) = ndata - 1
        ddymax = 0.0
        For i = 0 To pltvalue.nda(0)
            pltvalue.datax(0, i) = dt * CDbl(i)
            pltvalue.datay(0, i) = org_wave(i)
            If ddymax < Math.Abs(org_wave(i)) Then ddymax = Math.Abs(org_wave(i))
        Next i
        Call DEC_PLTDATA1(pltdata, ddymax, org_t_st, org_t_ed)
        Call SAKUZU(pltdata, pltvalue)
        fnameFIG = System.IO.Path.GetDirectoryName(Label2.Text) & "\fig-0-wave.png"
        PictureBox2.Image.Save(fnameFIG, System.Drawing.Imaging.ImageFormat.Png)
        '+++++++++++++++++++++++++++++++++++++++++++++
        '模擬波時刻歴画像出力
        '+++++++++++++++++++++++++++++++++++++++++++++
        pltvalue.ncase = 0
        ReDim pltvalue.nda(pltvalue.ncase)
        ReDim pltvalue.datax(pltvalue.ncase, ndata - 1)
        ReDim pltvalue.datay(pltvalue.ncase, ndata - 1)
        ReDim pltvalue.scom(pltvalue.ncase)
        pltvalue.scom(0) = "模擬波形"
        pltvalue.nda(0) = ndata - 1
        ddymax = 0.0
        For i = 0 To pltvalue.nda(0)
            pltvalue.datax(0, i) = dt * CDbl(i)
            pltvalue.datay(0, i) = ddy(i)
            If ddymax < Math.Abs(ddy(i)) Then ddymax = Math.Abs(ddy(i))
        Next i
        Call DEC_PLTDATA1(pltdata, ddymax, org_t_st, org_t_ed)
        Call SAKUZU(pltdata, pltvalue)
        fnameFIG = System.IO.Path.GetDirectoryName(Label2.Text) & "\fig-1-wave.png"
        PictureBox2.Image.Save(fnameFIG, System.Drawing.Imaging.ImageFormat.Png)
        '+++++++++++++++++++++++++++++++++++++++++++++
        '加速度応答スペクトル
        '+++++++++++++++++++++++++++++++++++++++++++++
        pltvalue.ncase = 2
        ReDim pltvalue.nda(pltvalue.ncase)
        ReDim pltvalue.datax(pltvalue.ncase, nfold - 1)
        ReDim pltvalue.datay(pltvalue.ncase, nfold - 1)
        ReDim pltvalue.scom(pltvalue.ncase)
        k = 0
        pltvalue.scom(k) = "原種波形"
        pltvalue.nda(k) = nfold - 1
        For i = 0 To pltvalue.nda(k)
            pltvalue.datax(k, i) = TK(i)
            pltvalue.datay(k, i) = org_res(i)
        Next i
        k = 1
        pltvalue.scom(k) = "目標SP"
        pltvalue.nda(k) = tsp_nd - 1
        For i = 0 To pltvalue.nda(k)
            pltvalue.datax(k, i) = tsp_sec(i)
            pltvalue.datay(k, i) = tsp_acc(i)
        Next i
        k = 2
        pltvalue.scom(k) = "模擬波形"
        pltvalue.nda(k) = nfold - 1
        For i = 0 To pltvalue.nda(k)
            pltvalue.datax(k, i) = TK(i)
            pltvalue.datay(k, i) = resacc(i)
        Next i
        Call DEC_PLTDATA2(pltdata)
        Call SAKUZU(pltdata, pltvalue)
        fnameFIG = System.IO.Path.GetDirectoryName(Label2.Text) & "\fig-2-sp.png"
        PictureBox2.Image.Save(fnameFIG, System.Drawing.Imaging.ImageFormat.Png)

    End Sub


    Private Sub FFT(ByVal nn As Integer, ByRef xr() As Double, ByRef xi() As Double)
        '**************************************************
        '高速フーリエ変換・逆変換
        '**************************************************
        'nn  ：データ数（2の累乗）
        'xr()：入出力データ実数部
        'xi()：入出力データ虚数部
        Dim g As Integer : Dim h As Integer : Dim i As Integer : Dim j As Integer : Dim k As Integer
        Dim l As Integer : Dim m As Integer : Dim n As Integer : Dim p As Integer : Dim q As Integer
        Dim a As Double : Dim b As Double : Dim xd As Double
        Dim s() As Double
        Dim c() As Double

        n = nn
        ReDim s(CInt(n / 2))
        ReDim c(CInt(n / 2))
        i = 0 : j = 0 : k = 0 : l = 0 : p = 0 : h = 0 : g = 0 : q = 0
        m = CInt(System.Math.Log(CDbl(n)) / System.Math.Log(2.0) + 1.0)
        a = 0.0 : b = Math.PI * 2.0 / CDbl(n)
        For i = 0 To CInt(n / 2)
            s(i) = System.Math.Sin(a) : c(i) = System.Math.Cos(a) : a = a + b
        Next i
        l = n : h = 1
        For g = 1 To m
            l = CInt(l / 2) : k = 0
            For q = 1 To h
                p = 0
                For i = k To l + k - 1
                    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 Then
                        xr(j) = a : xi(j) = b
                    Else
                        xr(j) = a * c(p) + b * s(p) : xi(j) = b * c(p) - a * s(p)
                    End If
                    p = p + h
                Next i
                k = k + l + l
            Next q
            h = h + h
        Next g
        j = CInt(n / 2)
        For i = 1 To n - 1
            k = n
            If j < i Then
                xd = xr(i) : xr(i) = xr(j) : xr(j) = xd
                xd = xi(i) : xi(i) = xi(j) : xi(j) = xd
            End If
            k = CInt(k / 2)
            Do While j >= k
                j = j - k : k = CInt(k / 2)
                If k = 0 Then Exit Do
            Loop
            j = j + k
        Next i

    End Sub


    Private Sub ERES(ByVal nn As Integer, ByVal dt As Double, ByVal h As Double, ByRef ddy() As Double, ByVal nt As Integer, _
                 ByRef tk() As Double, ByRef resacc() As Double, ByRef resvel() As Double, ByRef resdis() As Double)

        '地震応答スペクトル
        'nn      :加速度時刻歴総数
        'dt      :時間間隔（入力値）
        'h       :減衰定数（入力値）
        'ddy()   :加速度時刻歴（入力値）
        'nt      :計算周期総数(データ格納0〜nt-1)（入力値）
        'tk()    :計算周期（入力値）
        'resacc():加速度応答スペクトル（計算値）
        'resvel():速度応答スペクトル（計算値）
        'resdis():変位応答スペクトル（計算値）

        Dim i As Integer
        Dim m As Integer
        '応答最大値(gal,kine,cm)
        Dim accmax As Double
        Dim velmax As Double
        Dim dismax As Double

        Dim w As Double
        Dim w2 As Double
        Dim hw As Double
        Dim wd As Double
        Dim wdt As Double
        Dim e As Double
        Dim cwdt As Double
        Dim swdt As Double

        Dim ss As Double
        Dim cc As Double
        Dim s1 As Double
        Dim c1 As Double
        Dim s2 As Double
        Dim c2 As Double
        Dim s3 As Double
        Dim c3 As Double
        Dim A11 As Double
        Dim A12 As Double
        Dim A21 As Double
        Dim A22 As Double
        Dim B11 As Double
        Dim B12 As Double
        Dim B21 As Double
        Dim B22 As Double
        Dim dxf As Double
        Dim xf As Double
        Dim ddym As Double
        Dim ddyf As Double
        Dim x As Double
        Dim dx As Double
        Dim ddx As Double
        Dim pi As Double = Math.PI

        For i = 0 To nt - 1
            w = 2.0 * pi / tk(i)
            w2 = w * w
            hw = h * w
            wd = w * System.Math.Sqrt(1.0 - h * h)
            wdt = wd * dt
            e = System.Math.Exp(-hw * dt)
            cwdt = System.Math.Cos(wdt)
            swdt = System.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 * System.Math.Abs(ddy(0)) * dt
            velmax = System.Math.Abs(ddy(0)) * dt
            dismax = 0.0
            dxf = -ddy(0) * dt
            xf = 0.0

            For m = 1 To nn - 1
                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 <= System.Math.Abs(ddx)) Then accmax = System.Math.Abs(ddx)
                If (velmax <= System.Math.Abs(dx)) Then velmax = System.Math.Abs(dx)
                If (dismax <= System.Math.Abs(x)) Then dismax = System.Math.Abs(x)
                dxf = dx
                xf = x
            Next m
            resacc(i) = accmax
            resvel(i) = velmax
            resdis(i) = dismax
        Next i

    End Sub

    Private Sub IACC(ByVal nn As Integer, ByVal dt As Double, ByRef ddy() As Double, ByRef dy() As Double, ByRef y() As Double, _
                     ByRef ddymax As Double, ByRef dymax As Double, ByRef ymax As Double)

        '加速度時刻歴数値積分
        'nn:加速度時刻歴データ総数
        'dt    :時間間隔（入力値）
        'ddy() :加速度時刻歴（入力値）
        'dy()  :速度時刻歴（計算出力）
        'y()   :変位時刻歴（計算出力値）
        'ddymax:加速度最大値（計算出力値）
        'dymax :速度最大値（計算出力値）
        'ymax  :変位最大値（計算出力値）

        Dim i As Integer

        ddymax = 0.0
        dymax = 0.0
        ymax = 0.0

        dy(0) = 0.0
        y(0) = 0.0
        For i = 1 To nn - 1
            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 < System.Math.Abs(ddy(i))) Then ddymax = System.Math.Abs(ddy(i))
            If (dymax < System.Math.Abs(dy(i))) Then dymax = System.Math.Abs(dy(i))
            If (ymax < System.Math.Abs(y(i))) Then ymax = System.Math.Abs(y(i))
        Next i

    End Sub

    Private Sub CRAC(ByVal nn As Integer, ByVal dt As Double, ByRef ddy() As Double, ByRef dy() As Double, ByRef y() As Double, _
                     ByVal ddymax As Double, ByVal dymax As Double, ByVal ymax As Double)

        '加速度時刻歴基線補正
        'nn:加速度時刻歴データ総数
        'dt    :時間間隔（入力値）
        'ddy() :加速度時刻歴（入力値・書換出力値）
        'dy()  :速度時刻歴（入力値）
        'y()   :変位時刻歴（入力値）
        'ddymax:加速度最大値（入力値）
        'dymax :速度最大値（入力値）
        'ymax  :変位最大値（入力値）

        Dim i As Integer
        Dim tt As Double
        Dim t As Double
        Dim sum As Double
        Dim a1 As Double
        Dim a0 As Double
        Dim acmax As Double
        Dim coef As Double

        tt = CDbl(nn - 1) * dt
        t = 0.0
        For i = 0 To nn - 1
            y(i) = y(i) * (3.0 * tt - 2.0 * t) * t * t
            t = t + dt
        Next i
        sum = (y(0) + y(nn - 1)) / 2.0
        For i = 1 To nn - 2
            sum = sum + y(i)
        Next i
        sum = sum * dt
        a1 = 28.0 / 13.0 / tt / tt * (2.0 * dy(nn - 1) - 15.0 / System.Math.Pow(tt, 5.0) * sum)
        a0 = dy(nn - 1) / tt - a1 / 2.0 * tt
        t = 0.0
        acmax = 0.0
        For i = 0 To nn - 1
            ddy(i) = ddy(i) - a0 - a1 * t
            If (acmax < System.Math.Abs(ddy(i))) Then acmax = System.Math.Abs(ddy(i))
            t = t + dt
        Next i
        coef = ddymax / acmax
        For i = 0 To nn - 1
            ddy(i) = ddy(i) * coef
        Next i

    End Sub


    Private Sub SAKUZU(ByRef pltdata As DRAWPARA, ByRef pltvalue As DRAWVALUE)
        Dim i As Integer
        Dim j As Integer
        Dim pi As Double = Math.PI
        Dim thePictBox As PictureBox = PictureBox1
        Dim bmp As Bitmap
        Dim g As Graphics

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

        '対数値処理
        If pltdata.slogx = "L" Then
            For i = 0 To pltvalue.ncase
                For j = 0 To pltvalue.nda(i)
                    If pltvalue.datax(i, j) < System.Math.Pow(10.0, pltdata.xmin) Then pltvalue.datax(i, j) = System.Math.Pow(10.0, pltdata.xmin)
                    pltvalue.datax(i, j) = System.Math.Log10(pltvalue.datax(i, j))
                Next j
            Next i
        End If
        If pltdata.slogy = "L" Then
            For i = 0 To pltvalue.ncase
                For j = 0 To pltvalue.nda(i)
                    If pltvalue.datay(i, j) < System.Math.Pow(10.0, pltdata.ymin) Then pltvalue.datay(i, j) = System.Math.Pow(10.0, pltdata.ymin)
                    pltvalue.datay(i, j) = System.Math.Log10(pltvalue.datay(i, j))
                Next j
            Next i
        End If

        For i = 0 To 1
            Select Case i
                Case 0 '画面描画
                    pltdata.kpt = 1
                    thePictBox = PictureBox1
                    thePictBox.Visible = True
                Case 1 '画像書き出し
                    pltdata.kpt = 2
                    thePictBox = PictureBox2
                    thePictBox.Visible = False
            End Select
            thePictBox.Size = New Size(pltdata.kpt * pltdata.HSIZE, pltdata.kpt * pltdata.VSIZE)
            bmp = New Bitmap(thePictBox.Width, thePictBox.Height)
            thePictBox.Image = bmp
            g = Graphics.FromImage(thePictBox.Image)
            g.FillRectangle(Brushes.White, 0, 0, thePictBox.Width, thePictBox.Height)
            Call PLOT(g, pltdata, pltvalue)
            g.Dispose()
        Next i
        PictureBox1.Left = 0
        PictureBox1.Top = ToolStrip1.Height
        Me.ClientSize = New Size(PictureBox1.Width, PictureBox1.Height + ToolStrip1.Height)

    End Sub


    Private Sub PLOT(ByVal g As Graphics, ByVal pltdata As DRAWPARA, ByRef pltvalue As DRAWVALUE)

        Dim i As Integer
        Dim j As Integer
        Dim kxi As Integer : Dim kxf As Integer
        Dim kyi As Integer : Dim kyf As Integer
        Dim xx As Double : Dim yy As Double
        Dim kxx1 As Integer : Dim kyy1 As Integer
        Dim kxx2 As Integer : Dim kyy2 As Integer
        Dim str As String
        Dim f As New Font("ＭＳ ゴシック", pltdata.kpt * 10)
        Dim TextSize1 As New System.Drawing.SizeF
        Dim TextSize2 As New System.Drawing.SizeF
        Dim twmax As Integer
        Dim ds As Integer = 40

        Select Case pltvalue.ncase
            Case 0
                twmax = 0
            Case Else
                twmax = 0
                For i = 0 To pltvalue.ncase
                    TextSize1 = g.MeasureString(pltvalue.scom(i), f)
                    If twmax < CInt(TextSize1.Width) Then twmax = CInt(TextSize1.Width)
                Next i
                twmax = twmax + (ds + 20) * pltdata.kpt
        End Select
        kxi = pltdata.kpt * pltdata.kxi0
        kxf = pltdata.kpt * pltdata.kxf0 - twmax
        kyi = pltdata.kpt * pltdata.kyi0
        kyf = pltdata.kpt * pltdata.kyf0

        '座標軸
        Dim LPen As New System.Drawing.Pen(System.Drawing.Color.Black)
        LPen.DashStyle = Drawing2D.DashStyle.Dot
        If pltdata.slogx = "N" Then Call DRN_XJIKU(g, f, LPen, pltdata.kpt, pltdata.xmin, pltdata.xmax, pltdata.dx, kxi, kxf, kyi, kyf)
        If pltdata.slogy = "N" Then Call DRN_YJIKU(g, f, LPen, pltdata.kpt, pltdata.ymin, pltdata.ymax, pltdata.dy, kxi, kxf, kyi, kyf)
        If pltdata.slogx = "L" Then Call DRL_XJIKU(g, f, LPen, pltdata.kpt, pltdata.xmin, pltdata.xmax, kxi, kxf, kyi, kyf)
        If pltdata.slogy = "L" Then Call DRL_YJIKU(g, f, LPen, pltdata.kpt, pltdata.ymin, pltdata.ymax, 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)
        kxx1 = CInt(kxi - pltdata.kpt * 10 - TextSize1.Width - TextSize2.Height)
        kyy1 = CInt((kyi + kyf) / 2 + TextSize2.Width / 2)
        Call INC_STR(g, f, str, kxx1, kyy1, -90)
        'ｘ軸名描画
        str = pltdata.sxjiku
        TextSize2 = g.MeasureString(str, f)
        kxx1 = CInt((kxi + kxf) / 2 - TextSize2.Width / 2)
        kyy1 = CInt(kyf + pltdata.kpt * 10 + TextSize1.Height)
        Call INC_STR(g, f, str, kxx1, kyy1, 0)
        'データを線で連結
        LPen.Width = 1
        LPen.DashStyle = Drawing2D.DashStyle.Solid
        LPen.Color = Color.Blue
        For i = 0 To pltvalue.ncase
            If 1 <= pltvalue.ncase Then
                Select Case i
                    Case 0 : LPen.Color = Color.Lime : LPen.Width = 1
                    Case 1 : LPen.Color = Color.Red : LPen.Width = 5
                    Case 2 : LPen.Color = Color.Blue : LPen.Width = 2
                End Select
            End If
            xx = pltvalue.datax(i, 0)
            yy = pltvalue.datay(i, 0)
            kxx1 = kxi + CInt((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin))
            kyy1 = kyf - CInt((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin))
            For j = 1 To pltvalue.nda(i)
                xx = pltvalue.datax(i, j)
                yy = pltvalue.datay(i, j)
                kxx2 = kxi + CInt((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin))
                kyy2 = kyf - CInt((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin))
                If kxi <= kxx1 And kxx1 <= kxf Then
                    If kyi <= kyy1 And kyy1 <= kyf Then
                        If kxi <= kxx2 And kxx2 <= kxf Then
                            If kyi <= kyy2 And kyy2 <= kyf Then
                                g.DrawLine(LPen, kxx1, kyy1, kxx2, kyy2)
                            End If
                        End If
                    End If
                End If
                kxx1 = kxx2
                kyy1 = kyy2
            Next j
        Next i
        '凡例
        Select Case pltvalue.ncase
            Case 0
                str = pltvalue.scom(0)
                TextSize1 = g.MeasureString(str, f)
                g.DrawString(str, f, Brushes.Black, kxf - TextSize1.Width - pltdata.kpt * 2, kyi - TextSize1.Height - pltdata.kpt * 2)
            Case Else
                For i = 0 To pltvalue.ncase
                    str = pltvalue.scom(i)
                    TextSize1 = g.MeasureString(str, f)
                    g.DrawString(str, f, Brushes.Black, kxf + pltdata.kpt * (ds + 20), kyi + TextSize1.Height * i)
                    Select Case i
                        Case 0 : LPen.Color = Color.Lime : LPen.Width = 1
                        Case 1 : LPen.Color = Color.Red : LPen.Width = 5
                        Case 2 : LPen.Color = Color.Blue : LPen.Width = 2
                    End Select
                    kxx1 = kxf + pltdata.kpt * 10
                    kxx2 = kxx1 + ds * pltdata.kpt
                    kyy1 = kyi + CInt(TextSize1.Height * i) + CInt(TextSize1.Height / 2)
                    kyy2 = kyy1
                    g.DrawLine(LPen, kxx1, kyy1, kxx2, kyy2)
                Next i
                str = "PGA（gal）"
                TextSize1 = g.MeasureString(str, f)
                kyy1 = kyi + CInt(TextSize1.Height * 5)
                g.DrawString(str, f, Brushes.Black, kxf + pltdata.kpt * 10, kyy1)
                For i = 0 To pltvalue.ncase
                    str = (10 ^ pltvalue.datay(i, 0)).ToString("0") & "：" & pltvalue.scom(i)
                    TextSize1 = g.MeasureString(str, f)
                    g.DrawString(str, f, Brushes.Black, kxf + pltdata.kpt * 10, kyy1 + TextSize1.Height * (i + 1))
                Next i
        End Select
        f.Dispose()
        LPen.Dispose()

    End Sub


    Private Sub INC_STR(ByVal g As Graphics, ByVal f As System.Drawing.Font, ByVal str As String, _
                        ByVal kxx As Integer, ByVal kyy As Integer, ByVal ang As Single)
        '軸名描画
        g.ScaleTransform(1.0, 1.0)                '横・縦の表示比率を設定
        g.TranslateTransform(kxx, kyy)            '表示位置の設定（表示位置を原点とする座標移動）
        g.RotateTransform(ang)                    '表示角度を指定
        g.DrawString(str, f, Brushes.Black, 0, 0) '描画実行
        g.ResetTransform()                        '単位行列にリセット
    End Sub


    Private Sub DRN_XJIKU(ByVal g As Graphics, ByVal f As System.Drawing.Font, ByVal LPen As System.Drawing.Pen, _
                          ByVal kpt As Integer, _
                          ByVal xmin As Double, ByVal xmax As Double, ByVal dx As Double, _
                          ByVal kxi As Integer, ByVal kxf As Integer, _
                          ByVal kyi As Integer, ByVal kyf As Integer)
        Dim i As Integer
        Dim ix As Integer
        Dim xx As Double
        Dim wv As Double
        Dim str As String
        Dim kxx As Integer
        Dim TextSize1 As New System.Drawing.SizeF

        '普通ｘ軸描画
        ix = CInt((xmax - xmin) / dx)
        For i = 0 To ix
            xx = xmin + CDbl(i) * dx
            kxx = kxi + CInt((xx - xmin) * (kxf - kxi) / (xmax - xmin))
            wv = xmin + CDbl(i) * dx
            str = wv.ToString("0")
            If CInt(dx * 1000.0) Mod 10 <> 0 Then str = wv.ToString("0.000")
            If CInt(dx * 1000.0) Mod 10 = 0 And CInt(dx * 100.0) Mod 10 <> 0 Then str = wv.ToString("0.00")
            If (CInt(dx * 1000.0) Mod 10 = 0 And CInt(dx * 100.0) Mod 10 = 0) And CInt(dx * 10.0) Mod 10 <> 0 Then 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 + kpt * 3)
        Next i
    End Sub


    Private Sub DRN_YJIKU(ByVal g As Graphics, ByVal f As System.Drawing.Font, ByVal LPen As System.Drawing.Pen, _
                          ByVal kpt As Integer, _
                          ByVal ymin As Double, ByVal ymax As Double, ByVal dy As Double, _
                          ByVal kxi As Integer, ByVal kxf As Integer, _
                          ByVal kyi As Integer, ByVal kyf As Integer)
        Dim i As Integer
        Dim iy As Integer
        Dim yy As Double
        Dim wv As Double
        Dim str As String
        Dim kyy As Integer
        Dim TextSize1 As New System.Drawing.SizeF

        '普通ｙ軸描画
        iy = CInt((ymax - ymin) / dy)
        For i = 0 To iy
            yy = ymin + CDbl(i) * dy
            kyy = kyf - CInt((yy - ymin) * (kyf - kyi) / (ymax - ymin))
            wv = ymin + CDbl(i) * dy
            str = wv.ToString("0")
            If CInt(dy * 1000.0) Mod 10 <> 0 Then str = wv.ToString("0.000")
            If CInt(dy * 1000.0) Mod 10 = 0 And CInt(dy * 100.0) Mod 10 <> 0 Then str = wv.ToString("0.00")
            If (CInt(dy * 1000.0) Mod 10 = 0 And CInt(dy * 100.0) Mod 10 = 0) And CInt(dy * 10.0) Mod 10 <> 0 Then 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 - kpt * 2, kyy - TextSize1.Height / 2)
        Next i
    End Sub


    Private Sub DRL_XJIKU(ByVal g As Graphics, ByVal f As System.Drawing.Font, ByVal LPen As System.Drawing.Pen, _
                          ByVal kpt As Integer, _
                          ByVal xmin As Double, ByVal xmax As Double, _
                          ByVal kxi As Integer, ByVal kxf As Integer, _
                          ByVal kyi As Integer, ByVal kyf As Integer)
        Dim i As Integer
        Dim j As Integer
        Dim xx As Double
        Dim str As String
        Dim kxx As Integer
        Dim TextSize1 As New System.Drawing.SizeF

        '対数ｘ軸描画
        For i = CInt(xmin) To CInt(xmax) - 1
            For j = 1 To 9
                xx = System.Math.Log(CDbl(j) * System.Math.Pow(10.0, CDbl(i))) / System.Math.Log(10.0)
                kxx = kxi + CInt((xx - xmin) * (kxf - kxi) / (xmax - xmin))
                g.DrawLine(LPen, kxx, kyi, kxx, kyf)
            Next j
        Next i
        For i = CInt(xmin) To CInt(xmax)
            xx = CDbl(i)
            kxx = kxi + CInt((xx - xmin) * (kxf - kxi) / (xmax - xmin))
            str = CStr(System.Math.Pow(10.0, CDbl(i)))
            TextSize1 = g.MeasureString(str, f)
            g.DrawString(str, f, Brushes.Black, kxx - TextSize1.Width / 2, kyf + kpt * 3)
        Next i
    End Sub


    Private Sub DRL_YJIKU(ByVal g As Graphics, ByVal f As System.Drawing.Font, ByVal LPen As System.Drawing.Pen, _
                          ByVal kpt As Integer, _
                          ByVal ymin As Double, ByVal ymax As Double, _
                          ByVal kxi As Integer, ByVal kxf As Integer, _
                          ByVal kyi As Integer, ByVal kyf As Integer)
        Dim i As Integer
        Dim j As Integer
        Dim yy As Double
        Dim str As String
        Dim kyy As Integer
        Dim TextSize1 As New System.Drawing.SizeF

        '対数ｙ軸描画
        For i = CInt(ymin) To CInt(ymax) - 1
            For j = 1 To 9
                yy = System.Math.Log(CDbl(j) * System.Math.Pow(10.0, CDbl(i))) / System.Math.Log(10.0)
                kyy = kyf - CInt((yy - ymin) * (kyf - kyi) / (ymax - ymin))
                g.DrawLine(LPen, kxi, kyy, kxf, kyy)
            Next j
        Next i
        For i = CInt(ymin) To CInt(ymax)
            yy = CDbl(i)
            kyy = kyf - CInt((yy - ymin) * (kyf - kyi) / (ymax - ymin))
            str = CStr(System.Math.Pow(10.0, CDbl(i)))
            TextSize1 = g.MeasureString(str, f)
            g.DrawString(str, f, Brushes.Black, kxi - TextSize1.Width - kpt * 2, kyy - TextSize1.Height / 2)
        Next i
    End Sub


End Class