Option Explicit On
Option Strict On
Public Class Form1
    '***********************************************************
    '複素数計算のための構造体宣言とFunctionプロシージャ開始
    '***********************************************************
    Private Structure Complex
        Dim Real As Double
        Dim Imag As Double
    End Structure

    Private Function C_DOUBLE(ByVal a As Double, ByVal b As Double) As Complex
        '複素数実数部と虚数部への数値の代入
        C_DOUBLE.Real = a
        C_DOUBLE.Imag = b
    End Function

    Private Function C_PLUS(ByVal z1 As Complex, ByVal z2 As Complex) As Complex
        '加算
        C_PLUS.Real = z1.Real + z2.Real
        C_PLUS.Imag = z1.Imag + z2.Imag
    End Function

    Private Function C_MINUS(ByVal z1 As Complex, ByVal z2 As Complex) As Complex
        '減算
        C_MINUS.Real = z1.Real - z2.Real
        C_MINUS.Imag = z1.Imag - z2.Imag
    End Function

    Private Function C_TIMES(ByVal z1 As Complex, ByVal z2 As Complex) As Complex
        '乗算
        C_TIMES.Real = z1.Real * z2.Real - z1.Imag * z2.Imag
        C_TIMES.Imag = z1.Real * z2.Imag + z2.Real * z1.Imag
    End Function

    Private Function C_DIVIDED(ByVal z1 As Complex, ByVal z2 As Complex) As Complex
        '除算
        If z2.Real = 0.0 And z2.Imag = 0.0 Then
            MessageBox.Show("分母が0です！")
            Me.Close()
        Else
            C_DIVIDED.Real = (z1.Real * z2.Real + z1.Imag * z2.Imag) / (z2.Real ^ 2 + z2.Imag ^ 2)
            C_DIVIDED.Imag = (z2.Real * z1.Imag - z1.Real * z2.Imag) / (z2.Real ^ 2 + z2.Imag ^ 2)
        End If
    End Function

    Private Function D_ABS(ByVal z As Complex) As Double
        '絶対値
        D_ABS = Math.Sqrt(z.Real ^ 2 + z.Imag ^ 2)
    End Function

    Private Function C_CONJ(ByVal z As Complex) As Complex
        '共役複素数
        C_CONJ.Real = z.Real
        C_CONJ.Imag = -z.Imag
    End Function

    Private Function C_SQRT(ByVal z As Complex) As Complex
        '平方根
        If z.Real = 0.0 And z.Imag = 0.0 Then
            C_SQRT.Real = 0.0
            C_SQRT.Imag = 0.0
        Else
            Dim rr As Double
            Dim phi As Double
            rr = Math.Sqrt(z.Real ^ 2 + z.Imag ^ 2)
            phi = Math.Acos(z.Real / rr)
            If z.Imag < 0.0 Then phi = 2.0 * Math.PI - phi 'phiの範囲を0〜2πとする（acosの戻り値は0〜π）
            C_SQRT.Real = rr ^ (0.5) * Math.Cos(0.5 * phi)
            C_SQRT.Imag = rr ^ (0.5) * Math.Sin(0.5 * phi)
        End If
    End Function

    Private Function C_EXP(ByVal z As Complex) As Complex
        'EXP関数
        C_EXP.Real = Math.Exp(z.Real) * Math.Cos(z.Imag)
        C_EXP.Imag = Math.Exp(z.Real) * Math.Sin(z.Imag)
    End Function
    '***********************************************************
    '複素数計算のための構造体宣言とFunctionプロシージャ終了
    '***********************************************************

    Private Sub CFFT(ByVal n As Integer, ByRef x() As Complex, ByVal ind As Integer)
        '***********************************************************
        '複素数計算による高速フーリエ変換
        '***********************************************************
        Dim temp As Complex
        Dim theta As Complex
        Dim j As Integer
        Dim k As Integer
        Dim m As Integer
        Dim kmax As Integer
        Dim istep As Integer
        Dim pi As Double = Math.PI

        j = 1
        For i = 1 To n
            If i < j Then
                temp = x(j - 1)
                x(j - 1) = x(i - 1)
                x(i - 1) = temp
            End If
            m = CInt(n / 2)
            Do
                If j <= m Then Exit Do
                j = j - m
                m = CInt(m / 2)
            Loop While m >= 2
            j = j + m
        Next i
        kmax = 1
        Do
            istep = kmax * 2
            For k = 1 To kmax
                theta = C_DOUBLE(0.0, pi * CDbl(ind * (k - 1)) / CDbl(kmax))
                For i = k To n Step istep
                    j = i + kmax
                    temp = C_TIMES(x(j - 1), C_EXP(theta))
                    x(j - 1) = C_MINUS(x(i - 1), temp)
                    x(i - 1) = C_PLUS(x(i - 1), temp)
                Next i
            Next k
            kmax = istep
        Loop While kmax < n

    End Sub


    Private Sub CFESP(ByVal nl As Integer, ByRef Gth() As Double, ByRef Guw() As Double, _
                      ByRef GG() As Double, ByRef Galpha() As Double, ByRef Gbeta() As Double, _
                      ByVal lobj As Integer, ByVal lref As Integer, ByVal df As Double, _
                      ByVal nn As Integer, ByRef zz() As Complex, _
                      ByVal ind1 As Integer, ByVal ind2 As Integer)
        '***********************************************************
        '地盤の周波数応答関数（複素数計算）
        '***********************************************************
        'nl：基盤層を含む層の個数
        'Gth()：層厚(m)
        'Guw()：単位体積重量(tf/m3)
        'GG()：せん断弾性係数(tf/m3)
        'Galpha()：散乱減衰係数
        'Gbeta()：材料減衰係数
        'lobj：対象層の番号
        'lref：基準層の番号
        'df：周波数応答関数データの振動数間隔(Hz)
        'nn：周波数応答関数データの個数
        'zz()：周波数応答関数(無次元複素数)
        'ind1：0=加速度/加速度応答，1=ひずみ/加速度応答
        'ind2：0=基準層が内部，1=基準層が開放

        Dim i As Integer
        Dim k As Integer
        Dim gb() As Complex
        Dim p() As Complex
        Dim r() As Complex
        Dim a() As Complex
        Dim b() As Complex
        Dim ea As Complex
        Dim eb As Complex
        Dim ec As Complex
        Dim zobj As Complex
        Dim zref As Complex
        Dim w As Double '円振動数
        Dim dw As Double '円振動数増分
        Dim h As Double '減衰定数
        Dim cwork As Complex
        Dim cwork1 As Complex
        Dim cwork2 As Complex
        Dim pi As Double = Math.PI

        ReDim gb(nl - 1)
        ReDim p(nl - 1)
        ReDim r(nl - 2)
        ReDim a(nl - 1)
        ReDim b(nl - 1)

        w = 0.0
        dw = 2.0 * pi * df
        Select Case ind1
            Case 0 : zz(0) = C_DOUBLE(1.0, 0.0)
            Case 1 : zz(0) = C_DOUBLE(0.0, 0.0)
        End Select
        a(0) = C_DOUBLE(1.0, 0.0)
        b(0) = C_DOUBLE(1.0, 0.0)

        For k = 2 To nn
            w = w + dw
            For i = 1 To nl
                h = DAMPFAC(w, Galpha(i - 1), Gbeta(i - 1))
                gb(i - 1) = C_TIMES(C_DOUBLE(GG(i - 1), 0.0), C_DOUBLE(1.0, 2.0 * h))
                cwork = C_DIVIDED(C_DOUBLE(Guw(i - 1) / 9.8, 0.0), gb(i - 1))
                p(i - 1) = C_TIMES(C_DOUBLE(0.0, 1.0), C_SQRT(cwork))
            Next i
            For i = 1 To nl - 1
                r(i - 1) = C_TIMES(C_DIVIDED(gb(i - 1), gb(i)), C_DIVIDED(p(i - 1), p(i)))
            Next i
            For i = 1 To nl
                If i = lobj Then
                    Select Case ind1
                        Case 0
                            zobj = C_PLUS(a(i - 1), b(i - 1))
                        Case 1
                            cwork = C_DIVIDED(C_TIMES(p(i - 1), C_DOUBLE(w * Gth(i - 1), 0.0)), C_DOUBLE(2.0, 0.0))
                            ea = C_EXP(cwork)
                            eb = C_DIVIDED(C_DOUBLE(1.0, 0.0), ea)
                            cwork1 = C_MINUS(C_TIMES(a(i - 1), ea), C_TIMES(b(i - 1), eb))
                            cwork2 = C_DIVIDED(p(i - 1), C_DOUBLE(w, 0.0))
                            cwork = C_TIMES(C_TIMES(cwork1, cwork2), C_DOUBLE(0.01, 0.0))
                            zobj = C_TIMES(C_DOUBLE(-1.0, 0), cwork)
                    End Select
                End If
                If i = lref Then
                    Select ind2
                        Case 0 : zref = C_PLUS(a(i - 1), b(i - 1))
                        Case 1 : zref = C_TIMES(C_DOUBLE(2.0, 0.0), a(i - 1))
                    End Select
                End If
                If i < nl Then
                    cwork = C_TIMES(p(i - 1), C_DOUBLE(w * Gth(i - 1), 0.0))
                    ec = C_EXP(cwork)
                    ea = C_TIMES(a(i - 1), ec)
                    eb = C_DIVIDED(b(i - 1), ec)
                    cwork1 = C_PLUS(C_DOUBLE(1.0, 0.0), r(i - 1))
                    cwork1 = C_TIMES(cwork1, ea)
                    cwork2 = C_MINUS(C_DOUBLE(1.0, 0.0), r(i - 1))
                    cwork2 = C_TIMES(cwork2, eb)
                    cwork = C_PLUS(cwork1, cwork2)
                    a(i) = C_DIVIDED(cwork, C_DOUBLE(2.0, 0.0))
                    cwork1 = C_MINUS(C_DOUBLE(1.0, 0.0), r(i - 1))
                    cwork1 = C_TIMES(cwork1, ea)
                    cwork2 = C_PLUS(C_DOUBLE(1.0, 0.0), r(i - 1))
                    cwork2 = C_TIMES(cwork2, eb)
                    cwork = C_PLUS(cwork1, cwork2)
                    b(i) = C_DIVIDED(cwork, C_DOUBLE(2.0, 0.0))
                End If
            Next i
            zz(k - 1) = C_DIVIDED(zobj, zref)
        Next k

    End Sub

    Private Function DAMPFAC(ByVal w As Double, ByVal alpha As Double, ByVal beta As Double) As Double
        '減衰定数
        DAMPFAC = alpha / w + beta
    End Function

    Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click
        Dim i As Integer
        Dim dt As Double
        Dim ndata As Integer
        Dim nn As Integer
        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 fname1 As String = ""
        Dim fname2 As String = ""
        Dim ddy() As Double '基準層入力加速度
        Dim ddymax As Double '基準層入力加速度最大値
        Dim x() As Complex '複素フーリエ係数
        Dim acc() As Double '対象層加速度
        Dim accmax As Double '対象層加速度最大値

        Dim nl As Integer = 6
        Dim Gth() As Double = {5.0, 17.0, 13.0, 31.0, 4.0, 0.0}
        Dim Guw() As Double = {1.98, 1.95, 1.94, 1.94, 1.94, 1.83}
        Dim GG() As Double = {2900, 8000, 13400, 17800, 24200, 22800}
        Dim Galpha() As Double = {2, 2, 2, 2, 2, 2}
        Dim Gbeta() As Double = {0.02, 0.032, 0.024, 0.014, 0.01, 0.01}
        Dim lref As Integer = 6 '基準層（波形入力層）
        Dim lobj As Integer = 1 '対象層（波形出力層）
        Dim zz() As Complex '周波数応答関数（複素数）
        Dim df As Double
        Dim nfold As Integer
        Dim ind As Integer

        'ind = 11 '基準層：内部，対象層：内部
        'ind = 21 '基準層：内部，対象層：開放（対象層は基準層より下位）
        ind = 12 '基準層：解放，対象層：内部

        '地震加速度波形読込み
        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 ddy(ndata - 1)
        For i = 1 To ndata
            dat = sr.ReadLine : sbuf = dat.Split(delim) : ddy(i - 1) = CDbl(sbuf(0))
        Next i
        sr.Close()
        ddymax = 0.0
        For i = 1 To ndata
            If ddymax < Math.Abs(ddy(i - 1)) Then ddymax = Math.Abs(ddy(i - 1))
        Next i

        nn = 2
        Do
            nn = nn * 2
        Loop While nn < ndata * 1
        ReDim x(nn - 1)
        ReDim acc(nn - 1)
        nfold = CInt(nn / 2)
        ReDim zz(nfold - 1)

        For i = 1 To nn 'データ数を2の階乗にセット
            x(i - 1) = C_DOUBLE(0.0, 0.0)
        Next
        For i = 1 To ndata '加速度波形セット
            x(i - 1) = C_DOUBLE(ddy(i - 1), 0.0)
        Next i
        Call CFFT(nn, x, -1) 'フーリエ変換
        For i = 1 To nn
            x(i - 1) = C_DIVIDED(x(i - 1), C_DOUBLE(CDbl(nn), 0)) '出力値をデータ数で除す
        Next i
        '地盤の周波数応答関数計算
        df = 1.0 / CDbl(nn) / dt
        Select Case ind
            Case 11
                Call CFESP(nl, Gth, Guw, GG, Galpha, Gbeta, lobj, lref, df, nfold, zz, 0, 0)
            Case 12
                Call CFESP(nl, Gth, Guw, GG, Galpha, Gbeta, lobj, lref, df, nfold, zz, 0, 1)
            Case 21
                Call CFESP(nl, Gth, Guw, GG, Galpha, Gbeta, lref, lobj, df, nfold, zz, 0, 1)
                For k = 1 To nfold
                    zz(k - 1) = C_DIVIDED(C_DOUBLE(1.0, 0.0), zz(k - 1))
                Next k
        End Select
        '対象層加速度波形計算準備（フーリエ係数と地盤周波数応答関数の積）
        x(0) = C_TIMES(x(0), zz(0))
        For k = 2 To nfold - 1
            x(k - 1) = C_TIMES(x(k - 1), zz(k - 1))
            x(nn - k + 1) = C_CONJ(x(k - 1))
        Next k
        x(nfold - 1) = C_TIMES(x(nfold - 1), zz(nfold - 1))
        Call CFFT(nn, x, +1) 'フーリエ逆変換による対象層加速度波形算出
        accmax = 0.0
        For i = 1 To nn
            acc(i - 1) = x(i - 1).Real
            If accmax < Math.Abs(acc(i - 1)) Then accmax = Math.Abs(acc(i - 1))
        Next i

        '結果出力
        If SaveFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fname2 = SaveFileDialog1.FileName
        sw = New System.IO.StreamWriter(fname2, False, System.Text.Encoding.Default)
        dat = "--,基準層,対象層" : sw.WriteLine(dat)
        Select Case ind
            Case 11 : dat = "--,内部,内部" : sw.WriteLine(dat)
            Case 12 : dat = "--,解放,内部" : sw.WriteLine(dat)
            Case 21 : dat = "--,内部,解放" : sw.WriteLine(dat)
        End Select
        dat = "Layer No."
        dat = dat & "," & lref.ToString
        dat = dat & "," & lobj.ToString
        sw.WriteLine(dat)
        dat = "Max.gal"
        dat = dat & "," & ddymax.ToString("0.000")
        dat = dat & "," & accmax.ToString("0.000")
        sw.WriteLine(dat)
        dat = "dt"
        dat = dat & "," & dt.ToString("0.000")
        dat = dat & "," & dt.ToString("0.000")
        sw.WriteLine(dat)
        dat = "data_Nums."
        dat = dat & "," & ndata.ToString("0")
        dat = dat & "," & nn.ToString("0")
        sw.WriteLine(dat)
        dat = "i,ddy,acc" : sw.WriteLine(dat)
        For i = 1 To nn
            If i <= ndata Then
                dat = (i - 1).ToString
                dat = dat & "," & ddy(i - 1).ToString("0.000")
                dat = dat & "," & acc(i - 1).ToString("0.000")
                sw.WriteLine(dat)
            Else
                dat = (i - 1).ToString
                dat = dat & ",0," & acc(i - 1).ToString("0.000")
                sw.WriteLine(dat)
            End If
        Next i
        sw.Close()

    End Sub
End Class