Option Explicit On
Option Strict On
Public Class Form1
    Private NN As Integer = 30

    Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click

        main_part()

    End Sub

    Private Sub main_part()

        Dim i As Integer
        Dim k As Integer
        Dim ndata As Integer
        Dim nnn As Integer

        Dim pp As Double
        Dim ph As Double
        Dim pw As Double
        Dim r As Double
        Dim H As Double
        Dim z As Double
        Dim gamma As Double
        Dim kh As Double
        Dim lambda As Double
        Dim theta As Double
        Dim sum As Double
        Dim dz As Double
        Dim sw As System.IO.StreamWriter
        Dim fnameW As String = ""
        Dim fname As String = ""
        Dim dir As String = "c:\vb2008test\vbHDPsurge\data-file\"
        Dim dat As String = ""
        Dim xx(360) As Double
        Dim yy(360) As Double
        Dim px(360) As Double
        Dim py(360) As Double
        Dim ratio As Double = 2.0

        gamma = 1.0      '水の単位体積重量
        kh = 0.15        '水平震度
        r = 5.0          'サージタンク内半径
        H = 10.0         '水深
        theta = 0.0      '加速度の方向から半時計回りの角度

        FolderBrowserDialog1.Description = "データ出力フォルダを選択してください"
        If FolderBrowserDialog1.ShowDialog() = DialogResult.OK Then
            dir = FolderBrowserDialog1.SelectedPath
        End If

        For nnn = 0 To 3
            Select Case (nnn)
                Case 0 : H = 5.0 : fname = "hdp_H050.csv"
                Case 1 : H = 7.5 : fname = "hdp_H075.csv"
                Case 2 : H = 10.0 : fname = "hdp_H100.csv"
                Case 3 : H = 15.0 : fname = "hdp_H150.csv"
            End Select

            fnameW = dir & "\" & fname
            sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default)
            dat = "#z,pp,ph,pw"
            sw.WriteLine(dat)
            dz = 0.5
            z = -dz
            Do
                z = z + dz 'サージタンク底面からの距離
                sum = 0.0
                For i = 0 To NN
                    lambda = (2.0 * CDbl(i) + 1.0) / 2.0 * Math.PI
                    sum = sum + Math.Pow(-1.0, CDbl(i)) / lambda * FI(r / H, lambda) * Math.Cos(lambda * z / H)
                Next i
                '地震時動水圧強度(速度ポテンシャル)
                pp = gamma * kh * r * Math.Cos(theta) * sum
                '地震時動水圧強度(Housner)
                ph = Math.Sqrt(3.0) * gamma * kh * H * ((H - z) / H - 0.5 * (H - z) * (H - z) / H / H) * Math.Tanh(Math.Sqrt(3.0) * r / H) * Math.Cos(theta)
                '地震時動水圧強度(Westergaard)
                pw = 7.0 / 8.0 * gamma * kh * Math.Sqrt(H * (H - z))
                dat = z.ToString("0.000")
                dat = dat + "," + pp.ToString("0.000")
                dat = dat + "," + ph.ToString("0.000")
                dat = dat + "," + pw.ToString("0.000")
                sw.WriteLine(dat)
                If H - dz < z And z < H + dz Then Exit Do
            Loop While z < H
            sw.Close()
        Next nnn

        '動水圧水平分布
        H = 15.0
        fname = "dist_H150.csv"
        fnameW = dir & "\" & fname
        sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default)
        dat = "#x,y,dx,dy"
        sw.WriteLine(dat)
        z = 0.0 'サージタンク底面からの距離
        sum = 0.0
        For i = 0 To NN
            lambda = (2.0 * CDbl(i) + 1.0) / 2.0 * Math.PI
            sum = sum + Math.Pow(-1.0, CDbl(i)) / lambda * FI(r / H, lambda) * Math.Cos(lambda * z / H)
        Next i
        k = 0
        For i = 0 To 360 Step 5
            theta = CDbl(i) / 180.0 * Math.PI
            '地震時動水圧強度(速度ポテンシャル)
            pp = gamma * kh * r * Math.Cos(theta) * sum
            xx(k) = r * Math.Cos(theta)
            yy(k) = r * Math.Sin(theta)
            px(k) = ratio * pp * Math.Cos(theta)
            py(k) = ratio * pp * Math.Sin(theta)
            dat = xx(k).ToString("0.000")
            dat = dat + "," + yy(k).ToString("0.000")
            dat = dat + "," + px(k).ToString("0.000")
            dat = dat + "," + py(k).ToString("0.000")
            sw.WriteLine(dat)
            k = k + 1
        Next i
        sw.Close()
        ndata = k - 1

        '動水圧分布出力
        fname = "dtb_H150.csv"
        fnameW = dir & "\" & fname
        sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default)
        dat = "#x,y,dx,dy"
        sw.WriteLine(dat)
        For k = 0 To ndata
            dat = xx(k).ToString("0.000")
            dat = dat + "," + yy(k).ToString("0.000")
            dat = dat + "," + px(k).ToString("0.000")
            dat = dat + "," + py(k).ToString("0.000")
            sw.WriteLine(dat)
        Next k
        sw.Close()

        '動水圧成分出力
        fname = "pxy_H150.csv"
        fnameW = dir + "\" + fname
        sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default)
        'x方向（加速度方向）
        dat = "#x-direction" : sw.WriteLine(dat)
        dat = "#xx,yy,px,py" : sw.WriteLine(dat)
        For k = 0 To ndata
            If 0.0 <= xx(k) Then
                dat = r.ToString("0.000")
                dat = dat + "," + yy(k).ToString("0.000")
                dat = dat + "," + px(k).ToString("0.000")
                dat = dat + "," + (0.0).ToString("0.000")
                sw.WriteLine(dat)
            End If
            If xx(k) <= 0.0 Then
                dat = (-r).ToString("0.000")
                dat = dat + "," + yy(k).ToString("0.000")
                dat = dat + "," + px(k).ToString("0.000")
                dat = dat + "," + (0.0).ToString("0.000")
                sw.WriteLine(dat)
            End If
        Next k
        'y方向（加速度直交方向）
        dat = "#y-direction" : sw.WriteLine(dat)
        dat = "#xx,yy,px,py" : sw.WriteLine(dat)
        For k = 0 To ndata
            If 0.0 <= yy(k) Then
                dat = xx(k).ToString("0.000")
                dat = dat + "," + r.ToString("0.000")
                dat = dat + "," + (0.0).ToString("0.000")
                dat = dat + "," + py(k).ToString("0.000")
                sw.WriteLine(dat)
            End If
            If yy(k) <= 0.0 Then
                dat = xx(k).ToString("0.000")
                dat = dat + "," + (-r).ToString("0.000")
                dat = dat + "," + (0.0).ToString("0.000")
                dat = dat + "," + py(k).ToString("0.000")
                sw.WriteLine(dat)
            End If
        Next k
        sw.Close()

        MessageBox.Show("完了", "完了通知")

    End Sub

    Private Function FI(ByVal rbh As Double, ByVal lambda As Double) As Double
        Dim frbh As Double
        frbh = 2.0 / rbh * FMBF(1, lambda * rbh) / (lambda * FMBF(0, lambda * rbh) - FMBF(1, lambda * rbh) / rbh)
        Return frbh
    End Function

    Private Function FMBF(ByVal nu As Integer, ByVal z As Double) As Double
        '第1種変形Bessel関数
        Dim sum As Double
        Dim i As Integer
        sum = 0.0
        For i = 0 To NN
            sum = sum + 1.0 / GAMMAF(CDbl(i + 1)) / GAMMAF(CDbl(nu + i + 1)) * Math.Pow(0.5 * z, CDbl(2 * i))
        Next i
        Return Math.Pow(0.5 * z, CDbl(nu)) * sum
    End Function

    Private Function GAMMAF(ByVal z As Double) As Double
        'ガンマ関数
        Dim a As Double
        a = z / Math.E * Math.Sqrt(z * Math.Sinh(1.0 / z) + 1.0 / 810.0 / Math.Pow(z, 6.0))
        Return Math.Sqrt(2.0 * Math.PI / z) * Math.Pow(a, z)
    End Function

End Class