Option Explicit On
Option Strict On

Public Class Form1
    Private Structure SAKUZU '作図用データ構造体
        Dim sxjiku As String 'ｘ軸名
        Dim syjiku As String 'ｙ軸名
        Dim xmin As Double   'ｘ軸最小値
        Dim xmax As Double   'ｘ軸最大値
        Dim dx As Double     'ｘ軸増分
        Dim ymin As Double   'ｙ軸最小値
        Dim ymax As Double   'ｙ軸最大値
        Dim dy As Double     'ｙ軸増分
        Dim kxi0 As Integer  'ｘ軸最小値座標（px）
        Dim kxf0 As Integer  'ｘ軸最大値座標（px）
        Dim kyi0 As Integer  'ｙ軸最小値座標（px）
        Dim kyf0 As Integer  'ｙ軸最大値座標（px）
        Dim HSIZE As Integer '画像横寸法
        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）
    End Structure


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


    Private Sub main()
        'Fourier級数展開

        Dim scom As String
        Dim mm As Integer
        Dim ndata As Integer
        Dim datax() As Double
        Dim datay() As Double
        Dim i As Integer
        Dim k 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 fnameR As String = ""
        Dim fnameW As String = ""
        Dim fnameP As String
        Dim dt As Double
        Dim xr() As Double
        Dim xi() As Double
        Dim ak() As Double
        Dim bk() As Double
        Dim func() As Double

        Dim pltdata As SAKUZU

        'pltdata初期化
        pltdata.sxjiku = ""  'ｘ軸名
        pltdata.syjiku = ""  'ｙ軸名
        pltdata.xmin = 0.0   'ｘ軸最小値
        pltdata.xmax = 0.0   'ｘ軸最大値
        pltdata.dx = 0.0     'ｘ軸増分
        pltdata.ymin = 0.0   'ｙ軸最小値
        pltdata.ymax = 0.0   'ｙ軸最大値
        pltdata.dy = 0.0     'ｙ軸増分
        pltdata.kxi0 = 0     'ｘ軸最小値座標（px）
        pltdata.kxf0 = 0     'ｘ軸最大値座標（px）
        pltdata.kyi0 = 0     'ｙ軸最小値座標（px）
        pltdata.kyf0 = 0     'ｙ軸最大値座標（px）
        pltdata.HSIZE = 10000    '画像横寸法
        pltdata.DHL = 70      '画像左余白（縦軸描画スペース：px）
        pltdata.DHR = 20      '画像右余白（px）
        pltdata.VSIZE = 500    '画像縦サイズ（px）
        pltdata.DVL = 50      '画像下余白（横軸描画スペース：px）
        pltdata.DVU = 30      '画像上余白（px）

        '*******************************************************
        '入力ファイル名指定
        '*******************************************************
        OpenFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory()
        If OpenFileDialog1.ShowDialog() = DialogResult.OK Then fnameR = OpenFileDialog1.FileName
        If SaveFileDialog1.ShowDialog() = DialogResult.OK Then fnameW = SaveFileDialog1.FileName

        fnameP = System.IO.Path.GetDirectoryName(fnameW)
        dat = System.IO.Path.GetFileNameWithoutExtension(fnameW)
        fnameP = fnameP & "\fig_" & dat & ".png"

        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) : pltdata.sxjiku = sbuf(0) : pltdata.syjiku = sbuf(1)
        dat = sr.ReadLine() : sbuf = dat.Split(delim) : pltdata.HSIZE = CInt(sbuf(0)) : pltdata.DHL = CInt(sbuf(1)) : pltdata.DHR = CInt(sbuf(2))
        dat = sr.ReadLine() : sbuf = dat.Split(delim) : pltdata.VSIZE = CInt(sbuf(0)) : pltdata.DVL = CInt(sbuf(1)) : pltdata.DVU = CInt(sbuf(2))
        dat = sr.ReadLine() : sbuf = dat.Split(delim) : pltdata.xmin = CDbl(sbuf(0)) : pltdata.xmax = CDbl(sbuf(1)) : pltdata.dx = CDbl(sbuf(2))
        dat = sr.ReadLine() : sbuf = dat.Split(delim) : pltdata.ymin = CDbl(sbuf(0)) : pltdata.ymax = CDbl(sbuf(1)) : pltdata.dy = CDbl(sbuf(2))
        dat = sr.ReadLine() : sbuf = dat.Split(delim) : ndata = CInt(sbuf(0)) : mm = CInt(sbuf(1))
        nn = 2
        Do
            nn = nn * 2
        Loop While nn < ndata * 1 '必要に応じて後続の0を追加
        ReDim xr(nn - 1)
        ReDim xi(nn - 1)
        ReDim ak(CInt(nn / 2))
        ReDim bk(CInt(nn / 2))
        ReDim datax(ndata - 1)
        ReDim datay(ndata - 1)
        ReDim func(ndata - 1)

        For i = 1 To nn 'データ数を2の累乗にセット
            xr(i - 1) = 0.0 : xi(i - 1) = 0.0
        Next i
        For i = 1 To ndata
            dat = sr.ReadLine
            sbuf = dat.Split(delim)
            datax(i - 1) = CDbl(sbuf(0))
            datay(i - 1) = CDbl(sbuf(1))
        Next i
        sr.Close()

        dt = datax(1) - datax(0)
        For i = 1 To ndata
            xr(i - 1) = datay(i - 1) - datay(ndata - 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 = 0 To CInt(nn / 2)
            ak(i) = 2.0 * xr(i)
            bk(i) = -2.0 * xi(i)
        Next i
        For i = 0 To ndata - 1
            func(i) = 0.5 * ak(0)
            For k = 1 To mm
                func(i) = func(i) + ak(k) * Math.Cos(2.0 * Math.PI * k * i / nn / dt)
                func(i) = func(i) + bk(k) * Math.Sin(2.0 * Math.PI * k * i / nn / dt)
            Next k
            func(i) = func(i) + datay(ndata - 1)
        Next i

        sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default)
        sw.WriteLine(scom)
        sw.WriteLine("ndata," & ndata.ToString("0"))
        sw.WriteLine("x,y,func")
        For i = 0 To ndata - 1
            dat = datax(i).ToString("E")
            dat = dat & "," & datay(i).ToString("E")
            dat = dat & "," & func(i).ToString("E")
            sw.WriteLine(dat)
        Next i
        sw.Close()

        Call OUTFIG(fnameW, fnameP, pltdata)

        MessageBox.Show("計算完了")

    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 OUTFIG(ByVal fnameW As String, ByVal fnameP As String, ByRef pltdata As SAKUZU)
        Dim i As Integer
        Dim bmp As Bitmap
        Dim g As Graphics
        Dim kpt As Integer
        Dim thePicBox As PictureBox = PictureBox1

        Dim sr As System.IO.StreamReader
        Dim dat As String
        Dim sbuf() As String
        Dim delim() As Char = {","c}
        Dim nd As Integer
        Dim datax() As Double
        Dim datay() As Double
        Dim func() As Double
        Dim scom As String

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

        sr = New System.IO.StreamReader(fnameW, System.Text.Encoding.Default)
        dat = sr.ReadLine : sbuf = dat.Split(delim) : scom = sbuf(0)
        dat = sr.ReadLine : sbuf = dat.Split(delim) : nd = CInt(sbuf(1)) - 1
        dat = sr.ReadLine
        ReDim datax(nd)
        ReDim datay(nd)
        ReDim func(nd)
        For i = 0 To nd
            dat = sr.ReadLine : sbuf = dat.Split(delim)
            datax(i) = CDbl(sbuf(0))
            datay(i) = CDbl(sbuf(1))
            func(i) = CDbl(sbuf(2))
        Next i
        sr.Close()

        For i = 0 To 1
            Select Case i
                Case 0 '画面出力
                    kpt = 1
                    thePicBox = PictureBox1
                    thePicBox.Visible = True
                Case 1 '画像出力（ファイル：fname3，png形式）
                    kpt = 1
                    thePicBox = PictureBox2
                    thePicBox.Visible = False
            End Select
            thePicBox.Size = New Size(kpt * pltdata.HSIZE, 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)
            Call PLOT(g, kpt, pltdata, nd, datax, datay, func, scom)
            g.Dispose()
        Next i

        PictureBox1.Left = 0
        PictureBox1.Top = ToolStrip1.Height
        Me.ClientSize = New Size(PictureBox1.Width, PictureBox1.Height + ToolStrip1.Height)

        PictureBox2.Image.Save(fnameP, System.Drawing.Imaging.ImageFormat.Png)

    End Sub


    Private Sub PLOT(ByVal g As Graphics, ByVal kpt As Integer, ByVal pltdata As SAKUZU, _
                 ByVal nd As Integer, ByRef datax() As Double, ByRef datay() As Double, ByRef func() 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 kxx As Integer : Dim kyy As Integer
        Dim kxx1 As Integer : Dim kyy1 As Integer
        Dim kxx2 As Integer : Dim kyy2 As Integer
        Dim str As String
        Dim dc As Integer
        Dim f As New Font("ＭＳ ゴシック", kpt * 14)
        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 DRN_XJIKU(g, f, LPen, kpt, pltdata.xmin, pltdata.xmax, pltdata.dx, kxi, kxf, kyi, kyf)
        Call DRN_YJIKU(g, f, LPen, kpt, pltdata.ymin, pltdata.ymax, pltdata.dy, kxi, kxf, kyi, kyf)
        '枠線描画
        g.DrawRectangle(New Pen(Color.Black, 2), kxi, kyi, kxf - kxi, kyf - kyi)
        'コメント描画
        str = scom
        TextSize1 = g.MeasureString(str, f)
        g.DrawString(str, f, Brushes.Black, kxf - TextSize1.Width - kpt * 2, kyi - TextSize1.Height - 2 * kpt)
        'ｙ軸名描画
        str = pltdata.syjiku
        TextSize1 = g.MeasureString("-1.0", f)
        TextSize2 = g.MeasureString(str, f)
        kxx = CInt(kxi - kpt * 0 - TextSize1.Width - TextSize2.Height)
        kyy = CInt((kyi + kyf) / 2 + TextSize2.Width / 2)
        Call INC_STR(g, f, str, kxx, kyy, -90)
        'ｘ軸名描画
        str = pltdata.sxjiku
        TextSize2 = g.MeasureString(str, f)
        kxx = CInt((kxi + kxf) / 2 - TextSize2.Width / 2)
        kyy = CInt(kyf + kpt * 5 + TextSize1.Height)
        Call INC_STR(g, f, str, kxx, kyy, 0)
        'データプロット
        If kyf - kyi < kxf - kxi Then
            dc = CInt(kpt * (kyf - kyi) / 100)
        Else
            dc = CInt(kpt * (kxf - kxi) / 100)
        End If
        For i = 0 To nd
            xx = datax(i)
            yy = datay(i)
            kxx = kxi + CInt((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin)) - CInt(dc / 2)
            kyy = kyf - CInt((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin)) - CInt(dc / 2)
            g.DrawEllipse(New Pen(Color.Red, 1), kxx, kyy, dc, dc)
        Next i
        '回帰曲線
        xx = datax(0)
        yy = func(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 = func(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, 1), kxx1, kyy1, kxx2, kyy2)
            kxx1 = kxx2
            kyy1 = kyy2
        Next i
        '残差プロット
        If kyf - kyi < kxf - kxi Then
            dc = CInt(kpt * (kyf - kyi) / 100)
        Else
            dc = CInt(kpt * (kxf - kxi) / 100)
        End If
        For i = 0 To nd
            xx = datax(i)
            yy = datay(i) - func(i)
            kxx = kxi + CInt((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin)) - CInt(dc / 2)
            kyy = kyf - CInt((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin)) - CInt(dc / 2)
            g.DrawEllipse(New Pen(Color.Green, 1), kxx, kyy, dc, dc)
        Next i

        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 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



End Class