Option Explicit On
Option Strict On
Public Class Form1
    Private Structure SAKUZUL
        Dim sxjiku As String
        Dim syjiku As String
        Dim xmin As Double
        Dim xmax As Double
        Dim ymin As Double
        Dim ymax As Double
        Dim kxi0 As Integer
        Dim kxf0 As Integer
        Dim kyi0 As Integer
        Dim kyf0 As Integer
    End Structure


    Private Sub ToolStripButton1_Click_1(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click
        Call main()
    End Sub

    Private Sub main()
        Dim sr As System.IO.StreamReader
        Dim sw As System.IO.StreamWriter
        Dim dat As String
        Dim sbuf() As String
        Dim delim() As Char = {","c}
        Dim scom As String
        Dim nn As Integer
        Dim nt As Integer
        Dim dt As Double
        Dim h As Double
        Dim ddymax As Double
        Dim dymax As Double
        Dim ymax As Double
        Dim tk() As Double
        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 datax() As Double
        Dim datay() As Double
        Dim fnameR As String = ""
        Dim fnameW As String = ""
        Dim i As Integer

        h = 0.05 '減衰定数
        nt = 200 '周期分割数設定
        ReDim tk(nt - 1) : ReDim datax(nt - 1) : ReDim datay(nt - 1)

        '加速度時刻歴データ読み込み
        OpenFileDialog1.InitialDirectory = My.Computer.FileSystem.CurrentDirectory
        If OpenFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fnameR = OpenFileDialog1.FileName
        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) : dt = CDbl(sbuf(1))
        dat = sr.ReadLine() : sbuf = dat.Split(delim) : nn = CInt(sbuf(1))
        ReDim ddy(nn - 1) : ReDim dy(nn - 1) : ReDim y(nn - 1)
        ReDim resacc(nn - 1) : ReDim resvel(nn - 1) : ReDim resdis(nn - 1)
        For i = 1 To nn
            dat = sr.ReadLine() : sbuf = dat.Split(delim) : ddy(i - 1) = CDbl(sbuf(0))
        Next i
        sr.Close()

        '応答値算出
        For i = 0 To nt - 1
            tk(i) = System.Math.Pow(10.0, Math.Log10(2.0 * dt) + (Math.Log10(10.0) - Math.Log10(2.0 * dt)) / CDbl(nt) * CDbl(i))
        Next i
        Call IACC(nn, dt, ddy, dy, y, ddymax, dymax, ymax) '加速度時刻歴数値積分
        Call CRAC(nn, dt, ddy, dy, y, ddymax, dymax, ymax) '加速度時刻歴基線補正
        Call ERES(nn, dt, h, ddy, nt, tk, resacc, resvel, resdis) '応答スペクトル計算
        For i = 0 To nt - 1
            datax(i) = tk(i) : datay(i) = resacc(i)
        Next i

        '応答加速度書き込み
        '(作図時にデータの対数をとる前に普通軸データでの書き込みを先行して行う)
        If SaveFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fnameW = SaveFileDialog1.FileName
        sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default)
        sw.WriteLine(scom)
        sw.WriteLine("ndata," & CStr(nt))
        sw.WriteLine("period(sec),acc(gal),vel(kine),dis(cm)")
        For i = 0 To nt - 1
            sw.WriteLine(tk(i).ToString("0.00000") & "," & resacc(i).ToString("0.000") & "," & resvel(i).ToString("0.000") & "," & resdis(i).ToString("0.000"))
        Next
        sw.Close()

        '応答スペクトル（加速度）画面表示
        Dim kpt As Integer
        Dim HSIZE As Integer = 800
        Dim DHL As Integer = 80
        Dim DHR As Integer = 20
        Dim VSIZE As Integer = 600
        Dim DVL As Integer = 50
        Dim DVU As Integer = 20
        Dim nd As Integer
        Dim thePicBox As PictureBox = PictureBox1
        Dim bmp As Bitmap
        Dim g As Graphics
        Dim pltdata As SAKUZUL

        pltdata.sxjiku = "周期 (sec)"
        pltdata.syjiku = "加速度応答スペクトル (gal)"
        pltdata.xmin = -2.0
        pltdata.xmax = 1.0
        pltdata.ymin = 0.0
        pltdata.ymax = 4.0
        pltdata.kxi0 = DHL
        pltdata.kxf0 = HSIZE - DHR
        pltdata.kyi0 = DVU
        pltdata.kyf0 = VSIZE - DVL
        nd = nt - 1

        For i = 0 To nd
            datax(i) = Math.Log10(datax(i))
            datay(i) = Math.Log10(datay(i))
        Next i

        For i = 0 To 1
            Select Case i
                Case 0
                    kpt = 1
                    thePicBox = PictureBox1
                    thePicBox.Visible = True
                Case 1
                    kpt = 2
                    thePicBox = PictureBox2
                    thePicBox.Visible = False
            End Select
            thePicBox.Size = New Size(kpt * HSIZE, kpt * 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)
            Call PLOT(g, kpt, pltdata, nd, datax, datay, scom)
            g.Dispose()
        Next i
        '画面表示
        PictureBox1.Left = 0
        PictureBox1.Top = ToolStrip1.Height
        Me.ClientSize = New Size(PictureBox1.Width, PictureBox1.Height + ToolStrip1.Height)

        '画像保存
        fnameW = System.IO.Path.GetDirectoryName(fnameW) & "\" & _
          System.IO.Path.GetFileNameWithoutExtension(fnameW) & "-fig.png"
        PictureBox2.Image.Save(fnameW, System.Drawing.Imaging.ImageFormat.Png)

    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 PLOT(ByVal g As Graphics, ByVal kpt As Integer, ByVal pltdata As SAKUZUL, _
               ByVal nd As Integer, ByRef datax() As Double, ByRef datay() As Double, ByVal scom As String)

        Dim i 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("ＭＳ ゴシック", kpt * 12)
        Dim TextSize1 As New System.Drawing.SizeF
        Dim TextSize2 As New System.Drawing.SizeF

        kxi = kpt * pltdata.kxi0
        kxf = kpt * pltdata.kxf0
        kyi = kpt * pltdata.kyi0
        kyf = kpt * pltdata.kyf0

        '座標軸（対数ｘ軸，対数ｙ軸）
        Dim LPen As New System.Drawing.Pen(System.Drawing.Color.Black)
        LPen.DashStyle = Drawing2D.DashStyle.Dot
        Call DRL_XJIKU(g, f, LPen, kpt, pltdata.xmin, pltdata.xmax, kxi, kxf, kyi, kyf)
        Call DRL_YJIKU(g, f, LPen, 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 - 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 + kpt * 10 + TextSize1.Height)
        Call INC_STR(g, f, str, kxx1, kyy1, 0)
        'データを線で連結
        xx = datax(0)
        yy = datay(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 i = 1 To nd
            xx = datax(i)
            yy = datay(i)
            kxx2 = kxi + CInt((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin))
            kyy2 = kyf - CInt((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin))
            g.DrawLine(New Pen(Color.Blue, 2), kxx1, kyy1, kxx2, kyy2)
            kxx1 = kxx2
            kyy1 = kyy2
        Next i
        'コメント描画
        str = scom
        TextSize1 = g.MeasureString(str, f)
        g.DrawString(str, f, Brushes.Black, kxf - TextSize1.Width, kyi - TextSize1.Height - kpt * 1)
        f.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 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