Attribute VB_Name = "Module5"
'******************************************************************************
'*                                  マクロ                                    *
'******************************************************************************
'
'コンフィギュレーション識別番号：Mac-Ver18.0
'  配布用として不要なマクロ削除
'　2018/11/27　　原作者
'
' DataSmoothing                  ：データを平滑化する
' Demo_AttitudeIntegrations      ：姿勢パラメータ積分の演習（四元数、方向余弦マトリクス、オイラー角）
'       _Derivs                         微係数評価
' Demo_RandomNumber              ：正規分布の演習（正規分布をする２つの確率変数の和、差、積の分布）
' TableAndIntegTest              ：汎用テーブルと積分関数のテスト
'       _Derivs                         微係数評価
' TableTest                      ：汎用テーブルのテスト（マクロ自身で積分する例）


Option Explicit
Sub DataSmoothing()
'   DataSmoothing(Args(1), Args(2), Args(3), Args(4), Args(5))

'データを平滑化する

' Args(1) ：入力、独立変数データ（nR,1）
' Args(2) ：入力、従属変数データ（nR,nY）
' Args(3) ：入力、多項式近似次数（入力値）
' Args(4) ：入力、平滑化区間分割データ個数（１度に多項式近似するデータ個数）
' Args(5) ：入力、平滑化したデータの格納セル範囲（nR,nY）

    Dim L As Boolean        '論理型関数の戻り値
    Dim cSx As CellRange    '独立変数データのセル範囲（nR,1）
    Dim cSy As CellRange    '従属変数データのセル範囲（nR,nY）
    Dim cSz As CellRange    '平滑化したデータの格納セル範囲（nR,nY）
    Dim sx() As Double      '独立変数データ１次エリア（１次元）
    Dim sy() As Double      '従属変数データ１次エリア（２次元）
    Dim x() As Double       '独立変数データ多項式近似処理用エリア（１次元）
    Dim y() As Double       '従属変数データ多項式近似処理用エリア（２次元）
    Dim rS() As Double      '多項式近似残差
    Dim c() As Double       '多項式近似係数（２次元）
    Dim cF() As VectorND    '多項式近似係数（従属変数毎）
    Dim ix As Long          '独立変数データの列番号（従属変数列番号からの相対番号）
    Dim iy As Long          '従属変数データの列番号
    Dim iz As Long          '平滑化データの格納列番号
    Dim ir As Long          '最初のデータ行番号
    Dim nR As Long          'データ行数
    Dim nY As Integer       '従属変数の数
    Dim nS As Integer       '多項式近似スパン（１度に多項式近似するデータ個数）
    Dim nPi As Integer      '多項式近似次数（入力値）
    Dim nP As Integer       '多項式近似次数
    Dim nP1 As Integer      '多項式近似次数＋１
    Dim nRS As Integer      '多項式近似処理用データの外挿開始点
    Dim nSH As Integer      '多項式近似スパンの中間点
    Dim nC As Long          'シートからの読込バイアス行数
    Dim i As Long, j As Long, k As Long 'ループインディックス
    Dim n As Long, ip As Long  'ループインディックス
    Dim Yd As Double, rD As Double
    Dim aV As Double, dV As Double, Dvp As Double

    cSx = Args(1)                   '独立変数データのセル範囲
    cSy = Args(2)                   '従属変数データのセル範囲
    L = GetVariant(nPi, Args(3))    '平滑化の近似次数
    L = GetVariant(nS, Args(4))     '平滑化区間分割データ個数
    cSz = Args(5)                   '平滑化したデータの格納セル範囲

    iy = cSy.Column                 '従属変数データの先頭を列の基準とする
    nY = cSy.ColumnCount
    iz = cSz.Column - iy
    ix = cSx.Column - iy + 1
    ir = cSx.Row
    nR = cSx.RowCount

    'エリアの確保
    nSH = Int(nS / 2)
    nC = nS - nSH - 1
    nRS = nR - nC
    ReDim sx(nS)
    ReDim sy(nS, nY)
    ReDim x(nS)
    ReDim y(nS, nY)
    ReDim rS(nS)

    '多項式近似処理用１次エリアの初期設定
    Range("A1").Cells(ir, iy).Select
    For j = 1 To nS - nSH
        sx(j + nSH) = Selection.Cells(j, ix).Value
        For k = 1 To nY
            sy(j + nSH, k) = Selection.Cells(j, k).Value
        Next k
    Next j
    For j = 1 To nSH
        sx(j) = sx(nSH + 1) + (sx(nSH + 1) - Selection.Cells(2 * nSH + 2 - j, ix).Value)
        For k = 1 To nY
            sy(j, k) = sy(nSH + 1, k) + (sy(nSH + 1, k) - Selection.Cells(2 * nSH + 2 - j, k).Value)
        Next k
    Next j

    '平滑化処理
    For i = 1 To nR

        '１次エリアの設定
        If (i <= nRS) Then  '多項式近似スパンの最後にシートからデータを設定する
            sx(nS) = Selection.Cells(i + nC, ix).Value
            For k = 1 To nY
                sy(nS, k) = Selection.Cells(i + nC, k).Value
            Next k
        Else                '多項式近似スパン分のデータが無くなったので、データを外挿して追加する
            sx(nS) = sx(nS - 1) + (sx(nS - 1) - sx(nS - 2))
            For k = 1 To nY
                For j = 1 To nP1
                    cF(k).a(j) = c(j, k)
                Next j
                sy(nS, k) = PolyEx(cF(k).a, nP, sx(nS))
            Next k
        End If
        
        '１次エリアから多項式近似処理用エリアへデータを設定
        For j = 1 To nS
            x(j) = sx(j)
            For k = 1 To nY
                y(j, k) = sy(j, k)
            Next k
        Next j

        '残差の分散が最小となる次数を見つける
        Dvp = vbInf
        For ip = 1 To nPi

            '多項式近似処理用エリアの確保
            nP = ip
            nP1 = nP + 1
            ReDim c(nP1, nY)
            ReDim cF(nY)
            For k = 1 To nY
                ReDim cF(k).a(nP1)
            Next k

            '多項式近似処理
            L = LsqPol(x, y, nS, nY, c, nP)
            For k = 1 To nY
                For j = 1 To nP1
                    cF(k).a(j) = c(j, k)
                Next j
            Next k

            '残差の平均及び分散を求める
            aV = 0#
            dV = 0#
            For j = 1 To nS
                Yd = 0#
                rD = 0#
                For k = 1 To nY
                    rD = rD + PolyEx(cF(k).a, nP, x(j)) ^ 2
                    Yd = Yd + y(j, k) ^ 2
                Next k
                rS(j) = Sqr(rD) - Sqr(Yd)
                aV = aV + rS(j)
                dV = dV + rS(j) ^ 2
            Next j
            aV = aV / nS                    '残差の平均
            dV = dV / (nS - 1) - aV ^ 2     '残差の分散
            If (dV >= Dvp) Then Exit For    '分散が大きくなる方向に転じたらループを抜ける
            Dvp = dV
        Next ip

        '残差が分散の３倍以上のデータを棄却する
        dV = 3 * dV
        For j = 1 To nS
            If (rS(j) ^ 2 < dV) Then
                For k = 1 To nY
                    y(j, k) = PolyEx(cF(k).a, nP, x(j))
                Next k
            End If
        Next j

        '最適次数での近似係数を求める
        If (ip > nPi) Then nP = nP - 1
        nP1 = nP + 1
        ReDim c(nP1, nY)
        L = LsqPol(x, y, nS, nY, c, nP)

        '平滑化したデータをシートに設定する
        For k = 1 To nY
            For j = 1 To nP1
                cF(k).a(j) = c(j, k)
            Next j
            Selection.Cells(i, iz + k).Value = PolyEx(cF(k).a, nP, Selection.Cells(i, ix).Value)
        Next k

        '１次エリアをシフトする
        For j = 2 To nS
            sx(j - 1) = sx(j)
            For k = 1 To nY
                sy(j - 1, k) = sy(j, k)
            Next k
        Next j
    Next i

End Sub
Sub Demo_AttitudeIntegrations()
'   Demo_AttitudeIntegrations(Args(1), Args(2), Args(3))

'姿勢パラメータ積分の演習（レート積分、四元数、方向余弦マトリクス、オイラー角）

' Args(1)  ：入力、プログラムレートテーブルのセル範囲
' Args(2)  ：入力、積分間隔のセル範囲
' Args(3)  ：入力、積分する時刻のセル範囲（時刻の次のカラムから演算結果を格納する）

    Dim L As Boolean
    Dim cResult As CellRange                    '結果を格納するセル範囲
    Dim d() As Variant
    Dim dt As Double, TL As Double, Ta As Double
    Dim i As Integer, k As Integer

    L = GetTable(Table1, Args(1))               'プログラムレートテーブルの入力
    L = GetVariant(dt, Args(2))                 '積分間隔の入力

    cResult.Row = Args(3).Row                   '時刻セルの先頭行
    cResult.Column = Args(3).Column             '時刻セルの先頭列
    cResult.RowCount = 1                        '１行づつ処理する
    cResult.ColumnCount = Args(3).ColumnCount   'データ領域の列数

    L = GetArray(d, cResult)                    '初期値の取得
    INTG.nInt = Args(3).ColumnCount - 1

    L = IntIni(INTG.nInt)                       '積分領域の確保
    INTG.ID = 1                                 '微係数計算関数の識別
    INTG.t = d(1)                               '時刻
    INTG.dti = dt                               '通常の積分刻み
    
    For i = 1 To INTG.nInt                      '積分変数の初期値設定
        ORDN(i) = d(i + 1)
    Next i

    TL = Selection.Cells(Args(3).RowCount, 1).Value   'シートから計算終了時刻の取込み

    Do While INTG.t < TL                        '終了時刻まで繰り返す
        cResult.Row = cResult.Row + 1
        L = GetVariant(Ta, cResult)                 '計算時刻Taの取込み

        Do While (INTG.t < Ta)                      'Taまで積分する
            INTG.jIms = 0                               '不連続点での微係数を評価する
            L = Derivs
            INTG.tDis = TL                              '次の不連続点を設定する
            If (INTG.t < Table1.tDis And Table1.tDis < INTG.tDis) Then INTG.tDis = Table1.tDis
            If (INTG.t < Ta And Ta < INTG.tDis) Then INTG.tDis = Ta
    
            Call RkgInt                                 '不連続点(INTG.Tdis)まで積分する
        Loop

        d(1) = INTG.t                               '積分結果をシートに格納する
        For i = 1 To INTG.nInt
            d(i + 1) = ORDN(i)
        Next i
        L = PutArray(d, cResult)                '
    Loop

End Sub
Function Demo_AttitudeIntegrations_Derivs() As Boolean

'姿勢パラメータ積分の演習（四元数、方向余弦マトリクス、オイラー角）

    Dim wB As Vector3D, i As Integer, k As Integer
    Dim Qtn As Quaternion, Dcm As Matrix3x3, Ang As Vector3D
    Dim dQtn As Quaternion, dDcm As Matrix3x3, dAng As Vector3D

    'プログラムレート
    wB.x = TblInt(Table1, INTG.t)  'ロールレート
    wB.y = Table1.Val(2)           'ピッチレート
    wB.z = Table1.Val(3)           'ヨーレート
    DRVS(1) = wB.x
    DRVS(2) = wB.y
    DRVS(3) = wB.z

    '四元数
    Qtn.c = ORDN(4)
    Qtn.s.x = ORDN(5)
    Qtn.s.y = ORDN(6)
    Qtn.s.z = ORDN(7)
    dQtn = DrvQtn(Qtn, wB)
    DRVS(4) = dQtn.c
    DRVS(5) = dQtn.s.x
    DRVS(6) = dQtn.s.y
    DRVS(7) = dQtn.s.z

    '方向余弦マトリクス
    Dcm.x.x = ORDN(8)
    Dcm.x.y = ORDN(9)
    Dcm.x.z = ORDN(10)
    Dcm.y.x = ORDN(11)
    Dcm.y.y = ORDN(12)
    Dcm.y.z = ORDN(13)
    Dcm.x = VctNrm(Dcm.x)
    Dcm.y = VctNrm(Dcm.y)
    Dcm.z = VctCrp(Dcm.x, Dcm.y)
    ORDN(14) = Dcm.z.x
    ORDN(15) = Dcm.z.y
    ORDN(16) = Dcm.z.z
    dDcm = DrvDcm(Dcm, wB)
    DRVS(8) = dDcm.x.x
    DRVS(9) = dDcm.x.y
    DRVS(10) = dDcm.x.z
    DRVS(11) = dDcm.y.x
    DRVS(12) = dDcm.y.y
    DRVS(13) = dDcm.y.z

    'オイラー角（６種類）
    i = 16
    For k = 1 To 6
        Ang.x = ORDN(i + 1)
        Ang.y = ORDN(i + 2)
        Ang.z = ORDN(i + 3)
        dAng = DrvAng(Ang, wB, k)
        DRVS(i + 1) = dAng.x
        DRVS(i + 2) = dAng.y
        DRVS(i + 3) = dAng.z
        i = i + 3
    Next k

End Function
Sub Demo_RandomNumber()
'   Demo_RandomNumber(cRav1, cRsd1, cRav2, cRsd2, cCorr, cnR, cDst)

'正規分布の演習（正規分布をする２つの確率変数の和、差、積の分布）

' cRav1    ：入力、乱数１の平均のセル番地
' cRsd1    ：入力、乱数１の標準偏差のセル番地
' cRav2    ：入力、乱数２の平均のセル番地
' cRsd2    ：入力、乱数２の標準偏差のセル番地
' cCorr    ：入力、乱数１と乱数２の相関係数のセル番地
' cnR      ：入力、発生する乱数ペアの数のセル番地
' cDst     ：入力、平均、標準偏差、度数分布を格納する先頭行のセル範囲（nP行格納する）

    Dim L As Boolean
    Dim Rav1 As Double, Rsd1 As Double  '乱数１の平均､標準偏差
    Dim Rav2 As Double, Rsd2 As Double  '乱数２の平均､標準偏差
    Dim Corr As Double, nR As Long      '相関係数､発生する乱数ペアの数
    Dim cDst As CellRange               '平均、標準偏差、度数分布を格納するセル範囲
    Dim Pa(100) As Double               '２つの正規乱数の和の度数分布
    Dim Ps(100) As Double               '２つの正規乱数の和の度数分布
    Dim Pm(100) As Double               '２つの正規乱数の和の度数分布
    Dim i As Long, k As Long
    Dim Ru1 As Double, Ru2 As Double, Rup As Double
    Dim R1 As Double
    Dim R2 As Double
    Dim Cor2 As Double
    Dim rA As Double, RaAv As Double, RaSd As Double
    Dim rS As Double, RsAv As Double, RsSd As Double
    Dim Rm As Double, RmAv As Double, RmSd As Double
    Dim SaAv As Double, SaSd As Double
    Dim SsAv As Double, SsSd As Double
    Dim SmAv As Double, SmSd As Double
    Dim aScl As Double, sScl As Double, mScl As Double
    Dim nD As Long
    Dim nP As Long                      '度数分布を計数するエリアの数
    Dim nBas As Long                    '度数分布を計数するエリアの中間点

    L = GetVariant(Rav1, Args(1))       '引数の取り込み
    L = GetVariant(Rsd1, Args(2))
    L = GetVariant(Rav2, Args(3))
    L = GetVariant(Rsd2, Args(4))
    L = GetVariant(Corr, Args(5))
    L = GetVariant(nR, Args(6))
    cDst = Args(7)      '引数のセル番地の取り込み

    nP = 100            '度数分布を計数するエリアの数
    nBas = nP / 2
    Cor2 = Sqr(1 - Corr ^ 2)

    '２つの正規乱数の和
    RaAv = Rav1 + Rav2                          '作成する確率変数の平均（理論値）
    RaSd = Sqr(Rsd1 ^ 2 + 2# * Corr * Rsd1 * Rsd2 + Rsd2 ^ 2) '作成する確率変数の標準偏差（理論値）
    aScl = nP / (8# * RaSd)                     '度数分布計数用のスケーリング値
    SaAv = 0#                                   '作成した確率変数の平均
    SaSd = 0#                                   '作成した確率変数の標準偏差

    '２つの正規乱数の差
    RsAv = Rav1 - Rav2                          '作成する確率変数の平均（理論値）
    RsSd = Sqr(Rsd1 ^ 2 - 2# * Corr * Rsd1 * Rsd2 + Rsd2 ^ 2) '作成する確率変数の標準偏差（理論値）
    sScl = nP / (8# * RsSd)                     '度数分布計数用のスケーリング値
    SsAv = 0#                                   '作成した確率変数の平均
    SsSd = 0#                                   '作成した確率変数の標準偏差

    '２つの正規乱数の積
    RmAv = Corr * Rsd1 * Rsd2 + Rav1 * Rav2     '作成する確率変数の平均（理論値）
    RmSd = Sqr((Rsd1 * Rsd2) ^ 2 + 2 * (1 + Corr ^ 2) * Rsd1 * Rsd2) '作成する確率変数の標準偏差（理論値）
    mScl = nP / (8# * RmSd)                     '度数分布計数用のスケーリング値
    SmAv = 0#                                   '作成した確率変数の平均
    SmSd = 0#                                   '作成した確率変数の標準偏差

    For i = 1 To nR
        '相関係数を持つ２つの標準正規乱数
        Ru1 = rNd()
        Ru2 = rNd()
        Rup = Sqr(-2# * Log(Ru1))
        R1 = Rup * Sin(TOPI * Ru2)
        R2 = Corr * R1 + Cor2 * Rup * Cos(TOPI * Ru2)

        '変数（和、差、積）
        R1 = Rsd1 * R1 + Rav1
        R2 = Rsd2 * R2 + Rav2
        rA = R1 + R2
        rS = R1 - R2
        Rm = R1 * R2

        '変数の総和、二乗和
        SaAv = SaAv + rA
        SaSd = SaSd + rA ^ 2
        SsAv = SsAv + rS
        SsSd = SsSd + rS ^ 2
        SmAv = SmAv + Rm
        SmSd = SmSd + Rm ^ 2

        '度数分布の計数
        k = CInt(aScl * (rA - RaAv)) + nBas '和の分布
        If (k < 1) Then
            Pa(1) = Pa(1) + 1
        ElseIf (k > nP) Then
            Pa(nP) = Pa(nP) + 1
        Else
            Pa(k) = Pa(k) + 1
        End If
        k = CInt(sScl * (rS - RsAv)) + nBas '差の分布
        If (k < 1) Then
            Ps(1) = Ps(1) + 1
        ElseIf (k > nP) Then
            Ps(nP) = Ps(nP) + 1
        Else
            Ps(k) = Ps(k) + 1
        End If
        k = CInt(mScl * (Rm - RmAv)) + nBas '積の分布
        If (k < 1) Then
            Pm(1) = Pm(1) + 1
        ElseIf (k > nP) Then
            Pm(nP) = Pm(nP) + 1
        Else
            Pm(k) = Pm(k) + 1
        End If
    Next i
    SaAv = SaAv / nR                        '和の平均
    SaSd = Sqr(SaSd / (nR - 1) - SaAv ^ 2)  '和の標準偏差
    SsAv = SsAv / nR                        '差の平均
    SsSd = Sqr(SsSd / (nR - 1) - SsAv ^ 2)  '差の標準偏差
    SmAv = SmAv / nR                        '積の平均
    SmSd = Sqr(SmSd / (nR - 1) - SmAv ^ 2)  '積の標準偏差

    Range("A1").Cells(cDst.Row, cDst.Column).Select
    Selection.Cells(1, 1).Value = SaAv
    Selection.Cells(1, 2).Value = SaSd
    Selection.Cells(1, 3).Value = SsAv
    Selection.Cells(1, 4).Value = SsSd
    Selection.Cells(1, 5).Value = SmAv
    Selection.Cells(1, 6).Value = SmSd
    For k = 1 To nP
        Selection.Cells(k + 1, 1).Value = (k - nBas) / aScl + RaAv
        Selection.Cells(k + 1, 2).Value = Pa(k) / nR
        Selection.Cells(k + 1, 3).Value = (k - nBas) / sScl + RsAv
        Selection.Cells(k + 1, 4).Value = Ps(k) / nR
        Selection.Cells(k + 1, 5).Value = (k - nBas) / mScl + RmAv
        Selection.Cells(k + 1, 6).Value = Pm(k) / nR
    Next k

End Sub
Sub TableAndIntegTest()
'   TableAndIntegTest(Args(1), Args(2), Args(3))

' 汎用テーブルと積分関数のテスト

' Args(1)：入力、汎用関数テーブル
' Args(2)：入力、積分刻み、終了時間、出力時間間隔
' Args(3)：出力、出力領域（第1列は初期値が入っている）

    Dim Tin() As Variant, Din() As Variant, Dot(5) As Variant
    Dim TL As Double, Tp As Double, dTp As Double, L As Boolean
    Dim dtn As Double, dtd As Double
    Dim dtS As Double, dTdr As Double, dCor As Double
    Dim Rk1 As Double, Rk2 As Double, Rk3 As Double
    Dim i As Integer, k As Integer

    L = GetTable(Table1, Args(1))                           '汎用テーブル
    L = GetArray(Tin, Args(2))                              '積分刻み、終了時間、出力時間間隔
    L = GetArray(Din, Args(3))                              '初期値
    TL = Tin(2)
    dTp = Tin(3)
    
    '積分制御変数の初期化
    INTG.nInt = 2
    IntIni (INTG.nInt)
    INTG.t = Din(1)
    INTG.dti = Tin(1)
    INTG.dt = INTG.dti
    INTG.ID = 6                                             '微係数計算関数の識別
    
    '積分領域の初期化
    L = TableAndIntegTest_Derivs()
    ORDN(1) = Din(4)
    ORDN(2) = Din(5)
    
    '初期値の出力
    Tp = CDbl(CInt((INTG.t + vbTacr) / dTp + 1#)) * dTp
    Dot(1) = INTG.t
    Dot(2) = DRVS(1)
    Dot(3) = DRVS(2)
    Dot(4) = ORDN(1)
    Dot(5) = ORDN(2)
    L = PutArray(Dot, Args(3))
    
    Do While (INTG.t < TL)
        
        '次の不連続点の設定
        INTG.tDis = TL
        If (INTG.t < Table1.tDis And Table1.tDis < INTG.tDis) Then
            INTG.tDis = Table1.tDis
            INTG.jIms = 0
        End If
        If (INTG.t < Tp And Tp < INTG.tDis) Then INTG.tDis = Tp
        
        '不連続点(INTG.Tdis)まで積分
        Call RkgInt
                
        '計算結果の出力
        If (INTG.t - vbTacr < Tp And Tp < INTG.t + vbTacr) Then
            Tp = CDbl(CInt((Tp + vbTacr) / dTp + 1#)) * dTp
            Dot(1) = INTG.t
            Dot(2) = DRVS(1)
            Dot(3) = DRVS(2)
            Dot(4) = ORDN(1)
            Dot(5) = ORDN(2)
            Args(3).Row = Args(3).Row + 1
            L = PutArray(Dot, Args(3))
        End If
    Loop
    
End Sub
Function TableAndIntegTest_Derivs() As Boolean

    DRVS(1) = TblInt(Table1, INTG.t)
    DRVS(2) = Table1.Val(2)
    

End Function
Sub TableTest()
'   TableTest(Table1, argCell, outCell)

'汎用テーブルのテスト

' Table1  ：入力、テストするテーブル
' argCell ：入力、引数セル範囲
' outCell ：入力、出力セル範囲

    Dim L As Boolean, argCell As CellRange, outCell As CellRange
    Dim arg As Double, argp As Double, a As Double, F As Double
    Dim dltA As Double, dltT As Double, argT As Double, i As Integer
    ReDim Dot(14)
    
    Dim dtn As Double, dtd As Double      '時間間隔[s]
    Dim dtS As Double, dTdr As Double, dCor As Double   '積分で使用する変数
    Dim Rk1 As Double, Rk2 As Double, Rk3 As Double     'ルンゲ･クッタ･ギル積分定数の一時領域
    
    dltT = 2#

    L = GetTable(Table1, Args(1))           'テーブル

    argCell.Row = Args(2).Row               '引数セル位置の設定
    argCell.Column = Args(2).Column
    argCell.RowCount = 1
    argCell.ColumnCount = 1
    
    outCell.Row = Args(3).Row               '出力セル位置の設定
    outCell.Column = Args(3).Column
    outCell.RowCount = 1
    outCell.ColumnCount = Args(3).ColumnCount
    
    
    L = GetVariant(arg, argCell)            '最初の引数の入力
    a = arg
    
    '積分領域の初期化
    INTG.nInt = 2
    IntIni (INTG.nInt)
    INTG.t = 0#
    INTG.dti = 1#
    F = TblInt(Table1, a)
    DRVS(1) = 0#
    DRVS(2) = F
    ORDN(1) = a
    ORDN(2) = 0#
    
    '入力データがなくなるまで
    Do
        'データ入力
        argp = arg
        argCell.Row = argCell.Row + 1
        L = GetVariant(arg, argCell)        '次の引数の入力
        If (arg < vbTacr) Then Exit Do
        dltA = (arg - argp) / dltT          '引数の微係数はΔT間で一定とする
        argT = INTG.t + dltT
        
        '積分した引数が次の引数値に達するまで繰り返す
        Do While INTG.t < argT - vbTacr
            '不連続点処理
            INTG.jIms = 0
            a = ORDN(1)
            F = TblInt(Table1, a)
            DRVS(1) = dltA
            DRVS(2) = F
            
            '次の不連続点の設定
            INTG.tDis = INTG.t + dltT
            If (INTG.t < Table1.tDis - vbTacr And Table1.tDis < INTG.tDis) Then INTG.tDis = Table1.tDis
            
            '不連続点でのデータ出力
            Dot(1) = CDbl(INTG.jIms)
            Dot(2) = INTG.t
            Dot(3) = INTG.tDis
            Dot(4) = INTG.dt
            Dot(5) = a
            Dot(6) = F
            Dot(7) = ORDN(2)
            Dot(8) = CDbl(Table1.ix)
            Dot(9) = CDbl(Table1.ixp)
            Dot(10) = Table1.slpT
            Dot(11) = Table1.arg
            Dot(12) = Table1.aDis
            Dot(13) = Table1.slp
            Dot(14) = Table1.tDis
            L = PutArray(Dot, outCell)
            outCell.Row = outCell.Row + 1
            
            '次の不連続点まで積分
            Do While INTG.t < INTG.tDis - vbTacr
            
                dtn = CDbl(Int((INTG.t + vbTacr) / INTG.dti) + 1) * INTG.dti - INTG.t     '次の通常積分点までの時間
                dtd = INTG.tDis - INTG.t                                                    '不連続点までの時間
                If (vbTacr < dtd And dtd < dtn) Then INTG.dt = dtd Else INTG.dt = dtn
                
                For INTG.jIms = 1 To 4
                    Rk1 = RKG1(INTG.jIms)
                    Rk2 = RKG2(INTG.jIms)
                    Rk3 = RKG3(INTG.jIms)
                    dtS = RKG4(INTG.jIms) * INTG.dt
                    INTG.t = INTG.t + dtS
                    For i = 1 To INTG.nInt
                        dTdr = INTG.dt * DRVS(i)
                        dCor = Rk1 * (dTdr - Rk2 * TPI1(i))
                        ORDN(i) = ORDN(i) + dCor
                        TPI1(i) = TPI1(i) + 3# * dCor - Rk3 * dTdr
                        'a = ORDN(1)
                        'f = TblInt(table1, a)
                        'DRVS(2) = f
                    Next i
                    a = ORDN(1)
                    F = TblInt(Table1, a)
                    DRVS(2) = F
                    DRVS(1) = dltA
                    Dot(1) = CDbl(INTG.jIms)
                    Dot(2) = INTG.t
                    Dot(3) = INTG.tDis
                    Dot(4) = INTG.dt
                    Dot(5) = a
                    Dot(6) = F
                    Dot(7) = ORDN(2)
                    Dot(8) = CDbl(Table1.ix)
                    Dot(9) = CDbl(Table1.ixp)
                    Dot(10) = Table1.slpT
                    Dot(11) = Table1.arg
                    Dot(12) = Table1.aDis
                    Dot(13) = Table1.slp
                    Dot(14) = Table1.tDis
                    L = PutArray(Dot, outCell)
                    outCell.Row = outCell.Row + 1
                Next INTG.jIms
            Loop
        Loop
    Loop
End Sub
