円の内側の更に円の内側で円を転がす（式の一般化）

先日から、円の中で円を転がす多重スピログラフに挑戦している。
infoment.hatenablog.com
まとまった時間が取れない中で、さらに考えがまとまらず、結構な時間が
掛かってしまった。
今日も、前回の続きから。

横着して、今回はコメントを説明に代える。

Sub スピログラフ_2()

    Dim π As Double
        π = WorksheetFunction.Pi

    Dim r(1 To 4) As Long
        r(1) = 11
        r(2) = 5
        r(3) = 2
        r(4) = 2
    
    ' 円の分割数。
    ' 36であれば、36点で円を描画する。点から点の回転角は360/36=10°となる。
    Dim division_count As Long
        division_count = 36
    
    ' ２個めの円の回転回数。
    Dim RoundNumber As Long
    
    ' division_countを何回繰り返すか。
    Dim AllDivisionCount As Long
    Dim i As Long
    Dim j As Long
    Dim k As Long
        ' ２個めの円の回転回数。１個目と２個目の半径の最小公倍数から、
        ' ２個めの円が何回転すると元の位置に戻るかを求めている。
        RoundNumber = WorksheetFunction.Lcm(r(1), r(2)) / r(1)
        
        ' division_count、つまり描画点から次の描画点までの角度を何度
        ' 繰り返せば元の位置に戻るかを求めている。
        ' ※最後に１加えるのは、最初と最後の点を重ねて円が切れないようにするため。
        AllDivisionCount = division_count * RoundNumber + 1
            
    ' ｎ番目の円の中心座標。
    Dim ox() As Double: ReDim ox(1 To UBound(r), 1 To AllDivisionCount)
    Dim oy() As Double: ReDim oy(1 To UBound(r), 1 To AllDivisionCount)
    
        ' 一番大きな円は動かないので、中心は常に原点。
        For j = 1 To AllDivisionCount
            ox(1, j) = 0
            oy(1, j) = 0
        Next

    ' 表の最大列数。
    Dim cMax As Long
        cMax = UBound(r) + 1

    ' グラフリセット。
        Sheet3.Cells.ClearContents

    ' 各円を分割する角度（rad）。
    Dim dθ() As Double
    ReDim dθ(1 To UBound(r) - 1)
        For i = 1 To UBound(dθ)
            Select Case i
                Case 1
                    dθ(1) = 2 * π / division_count
                Case Else
                    dθ(i) = -r(i - 1) / r(i) * dθ(i - 1)
            End Select
        Next

    ' 全描画に必要なグラフの行数。
    ' ※円の描画に必要な点　division_count+1 ←始点と終点を重ねないと円が閉じない。
    ' ※動点Ｐの描画に必要な点　回転数の掛け算+1
    Dim rMax As Long
        rMax = (division_count + 1) * (UBound(r) - 1) + AllDivisionCount

    ' 各座標を格納するための配列。
    ' 一次元を０始まりにするのは、ラベル名を格納するため。
    Dim c As Long
    Dim arr() As Variant
    ReDim arr(rMax + 1, 1 To cMax + 1)
        arr(0, 1) = "X"
        For c = 2 To cMax + 1
            arr(0, c) = "Y" & c - 1
        Next

    Dim myChart As ChartObject
    Set myChart = ActiveSheet.ChartObjects(1)
        myChart.Chart.SetSourceData Source:=Sheet3.Cells(1, 1).Resize(rMax + 2, cMax + 1)


    ' 各座標。
    Dim x() As Double: ReDim x(1 To UBound(r), 1 To division_count + 1)
    Dim y() As Double: ReDim y(1 To UBound(r), 1 To division_count + 1)

    ' 以降、円を描画するための変数はｊとする。
    ' 以降、AllDivisionCount番目の中心座標を求めるための変数はkとする。
    
    ' 一番大きな円を描画するための、変化しない座標。
        For j = 1 To division_count + 1
            ' １個目の円の円周部分を描画する座標。
            x(1, j) = ox(1, j) + r(1) * Cos(dθ(1) * (j - 1))
            y(1, j) = oy(1, j) + r(1) * Sin(dθ(1) * (j - 1))
            arr(j, 1) = x(1, j)
            arr(j, 2) = y(1, j)
        Next
    
    ' 一番大きな円の描画に合わせて、二つ目以降の円の中心座標をAllDivisionCount分求める。
    ' ※下記の計算で、上記で計算したｘｙ座標を使用しないのは、上記の計算結果が１周分しかないため。
    ' 　途中で計算結果が無くなってしまうので（=division_count+2個め以降がない）、改めて計算している。
    Dim n As Long
        For k = 1 To AllDivisionCount
            For n = 2 To UBound(r)
                ox(n, k) = ox(n - 1, k) + r(n - 1) * Cos(dθ(n - 1) * (k - 1)) - r(n) * Cos(dθ(n - 1) * (k - 1))
                oy(n, k) = oy(n - 1, k) + r(n - 1) * Sin(dθ(n - 1) * (k - 1)) - r(n) * Sin(dθ(n - 1) * (k - 1))
            Next
        Next
    
    ' 以降、動く円の描画。
        For k = 1 To AllDivisionCount
            ' n個目の円の円周を描画する座標。
            For n = 2 To UBound(r) - 1
                For j = 1 To division_count + 1                    '
                    arr(j + (n - 1) * (division_count + 1), 1) = ox(n, k) + r(n) * Cos(dθ(n) * (j - 1))
                    arr(j + (n - 1) * (division_count + 1), n + 1) = oy(n, k) + r(n) * Sin(dθ(n) * (j - 1))
                Next
            Next
                
            ' 最後は円を描かず、その中心座標のみで描画（軌跡を描く）。
            arr(k + (n - 1) * (division_count + 1), 1) = ox(n, k)
            arr(k + (n - 1) * (division_count + 1), n + 1) = oy(n, k)
            
            ' 描画と一緒にポイントを動かす。
            arr(k + (n - 1) * (division_count + 1) + 1, 1) = ox(n, k)
            arr(k + (n - 1) * (division_count + 1), n + 2) = vbNullString
            arr(k + (n - 1) * (division_count + 1) + 1, n + 2) = oy(n, k)
            
    
            Application.Wait [Now() + "00:00:00.03"]
            DoEvents
            
            Sheet3.Cells(1, 1).Resize(rMax + 2, cMax + 1) = arr
        
        Next
            
End Sub

実行してみた。こんな感じだ。

次回に続きます。

参考まで。