からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

【R】黄金角の定義の可視化

はじめに

 円関数(三角関数)の定義や公式をR言語で可視化して理解しようシリーズの補足シリーズ(黄金比編)です。

 この記事では、黄金角の定義を確認します。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

黄金角の定義の可視化

 黄金比(golden ratio)を用いて定義される黄金角(golden angle)の定義をグラフで確認します。

 利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)
library(patchwork)

 この記事では、基本的に パッケージ名::関数名() の記法を使うので、パッケージの読み込みは不要です。ただし、作図コードについてはパッケージ名を省略するので、 ggplot2 を読み込む必要があります。
 また、ネイティブパイプ演算子 |> を使います。magrittr パッケージのパイプ演算子 %>% に置き換えられますが、その場合は magrittr を読み込む必要があります。

黄金角の定義

 黄金角の定義を数式とグラフで確認します。
 円周の座標や度数法と弧度法の角度については「【R】円周の作図 - からっぽのしょこ」を参照してください。

 黄金角は、円周を黄金比で分割した2つの円弧の内、小さい方の弧の中心角として定義されます。

・作図コード(クリックで展開)

 黄金数を作成します。

# 黄金数を計算
phi <- (1 + sqrt(5)) * 0.5

# 黄金数の逆数を計算
recip_phi <- 1 / phi
phi; recip_phi
[1] 1.618034
[1] 0.618034

 黄金数(黄金比の値)  \varphi = \frac{1 + \sqrt{5}}{2} と黄金数の逆数  \frac{1}{\varphi} を計算します。

 半径を指定して、円周と円弧の長さを作成します。

# 半径を指定
r <- 2

# 円周の長さを計算
ab <- 2 * pi * r

# 円弧の長さを計算
a <- ab * recip_phi
b <- ab * (1 - recip_phi)
ab; a; b
[1] 12.56637
[1] 7.766444
[1] 4.799926

 半径  r を指定して、円周  a + b = 2 \pi r と円孤  a = (a + b) \frac{1}{\varphi} b = (a + b) (1 - \frac{1}{\varphi}) を計算します。

 黄金角を作成します。

# 中心角を計算
#alpha <- a / r
#beta  <- b / r
alpha <- 2*pi * recip_phi
beta  <- 2*pi * (1 - recip_phi)
alpha; beta
[1] 3.883222
[1] 2.399963

 円弧  a, b それぞれの中心角  \alpha = \frac{a}{r} = 2 \pi \frac{1}{\varphi} \beta = \frac{b}{r} = 2 \pi (1 - \frac{1}{\varphi}) を計算します。

 円周の描画用のデータフレームを作成します。

# 円周の座標を作成
circle_df <- tibble::tibble(
  theta = seq(from = 0, to = 2*pi, length.out = 901), # ラジアン
  x     = r * cos(theta), 
  y     = r * sin(theta), 
  arc_label = dplyr::if_else(
    theta <= alpha, true = "a", false = "b"
  ) # 円弧ラベル
)
circle_df
# A tibble: 901 × 4
     theta     x      y arc_label
     <dbl> <dbl>  <dbl> <chr>    
 1 0        2    0      a        
 2 0.00698  2.00 0.0140 a        
 3 0.0140   2.00 0.0279 a        
 4 0.0209   2.00 0.0419 a        
 5 0.0279   2.00 0.0558 a        
 6 0.0349   2.00 0.0698 a        
 7 0.0419   2.00 0.0838 a        
 8 0.0489   2.00 0.0977 a        
 9 0.0559   2.00 0.112  a        
10 0.0628   2.00 0.126  a        
# ℹ 891 more rows

  0 \leq \theta \leq 2 \pi の範囲のラジアンを作成して、半径が  r の円周のx軸の値  x = r \cos \theta、y軸の値  y = r \sin \theta を計算します。
 弧  a の範囲  0 \leq \theta \leq \alpha と弧  b の範囲  \alpha \lt \theta \leq 2 \pi それぞれにラベル付けします。

 中心角の描画用のデータフレームを作成します。

# 角度線の座標を作成
radius_df <- tibble::tibble(
  theta = c(0, alpha), 
  x     = r * cos(theta), 
  y     = r * sin(theta), 
  edge_label = c("r", NA)
)
radius_df
# A tibble: 2 × 4
  theta     x     y edge_label
  <dbl> <dbl> <dbl> <chr>     
1  0     2     0    r         
2  3.88 -1.47 -1.35 <NA>    

 ラジアン  \theta = 0, \alpha における円周上の点の座標を計算します。

 角マークの描画用のデータフレームを作成します。

# 角マークの座標を作成
d_a <- 0.15
d_b <- 0.2
angle_arc_df <- dplyr::bind_rows(
  tibble::tibble(
    theta = seq(from = 0, to = alpha, length.out = 100), 
    x     = d_a * cos(theta), 
    y     = d_a * sin(theta), 
    arc_label = "a"
  ), 
  tibble::tibble(
    theta = seq(from = 0, to = -beta, length.out = 100), 
    x     = d_b * cos(theta), 
    y     = d_b * sin(theta), 
    arc_label = "b"
  )
)
angle_arc_df
# A tibble: 200 × 4
    theta     x       y arc_label
    <dbl> <dbl>   <dbl> <chr>    
 1 0      0.15  0       a        
 2 0.0392 0.150 0.00588 a        
 3 0.0784 0.150 0.0118  a        
 4 0.118  0.149 0.0176  a        
 5 0.157  0.148 0.0234  a        
 6 0.196  0.147 0.0292  a        
 7 0.235  0.146 0.0350  a        
 8 0.275  0.144 0.0407  a        
 9 0.314  0.143 0.0463  a        
10 0.353  0.141 0.0519  a        
# ℹ 190 more rows

 弧  a の範囲  0 \leq \theta \leq \alpha と弧  b の範囲  \alpha \lt \theta \leq 2 \pi のラジアンを用いて、それぞれの中心からの位置(半径)を d_a, d_b として、円周の座標と同様に計算します。弧  b のラジアンの範囲は、負の角度を用いて  -\beta \leq \theta \leq 0 でも計算できます。

 弧ラベルと角ラベルの描画用のデータフレームを作成します。

# 角ラベルの座標を作成
angle_label_df <- tibble::tibble(
  theta = 0.5 * c(alpha, -beta), 
  x     = cos(theta), 
  y     = sin(theta), 
  arc_label   = c("a", "b"), 
  angle_label = c("alpha", "beta")
)
angle_label_df
# A tibble: 2 × 5
  theta      x      y arc_label angle_label
  <dbl>  <dbl>  <dbl> <chr>     <chr>      
1  1.94 -0.362  0.932 a         alpha      
2 -1.20  0.362 -0.932 b         beta       

  a, b の弧ラベルと  \alpha, \beta の角ラベルをそれぞれの円弧の中点に表示することにします。
 ラジアン  \theta = \frac{\alpha}{2}, \alpha + \frac{\beta}{2} を用いて、半径が1の円弧上の点の座標を計算します。弧  b のラジアンは、負の角度を用いて  -\frac{\beta}{2} でも計算できます。

 角度目盛の描画用のデータフレームを作成します。

# 半円における目盛数(分母の値)を指定
denom <- 6

# 角度目盛線の座標を作成
angle_axis_df <- tibble::tibble(
  i     = 0:(2*denom-1),  # 目盛位置番号(分子の値)
  theta = i / denom * pi, # 目盛値
  x     = r * cos(theta), 
  y     = r * sin(theta), 
  rad_label = paste0("frac(", i, ", ", denom, ") ~ pi"), # 角度ラベル
  h = 1 - (x/r * 0.5 + 0.5), # ラベル位置の調整値
  v = 1 - (y/r * 0.5 + 0.5)  # ラベル位置の調整値
)
angle_axis_df
# A tibble: 12 × 7
       i theta         x         y rad_label             h      v
   <int> <dbl>     <dbl>     <dbl> <chr>             <dbl>  <dbl>
 1     0 0      2   e+ 0  0        frac(0, 6) ~ pi  0      0.5   
 2     1 0.524  1.73e+ 0  1   e+ 0 frac(1, 6) ~ pi  0.0670 0.25  
 3     2 1.05   1   e+ 0  1.73e+ 0 frac(2, 6) ~ pi  0.25   0.0670
 4     3 1.57   1.22e-16  2   e+ 0 frac(3, 6) ~ pi  0.5    0     
 5     4 2.09  -1   e+ 0  1.73e+ 0 frac(4, 6) ~ pi  0.75   0.0670
 6     5 2.62  -1.73e+ 0  1   e+ 0 frac(5, 6) ~ pi  0.933  0.25  
 7     6 3.14  -2   e+ 0  2.45e-16 frac(6, 6) ~ pi  1      0.5   
 8     7 3.67  -1.73e+ 0 -1   e+ 0 frac(7, 6) ~ pi  0.933  0.75  
 9     8 4.19  -1.00e+ 0 -1.73e+ 0 frac(8, 6) ~ pi  0.75   0.933 
10     9 4.71  -3.67e-16 -2   e+ 0 frac(9, 6) ~ pi  0.5    1     
11    10 5.24   1   e+ 0 -1.73e+ 0 frac(10, 6) ~ pi 0.25   0.933 
12    11 5.76   1.73e+ 0 -1.00e+ 0 frac(11, 6) ~ pi 0.0670 0.75  

 「円周の作図」を参照してください。

 円弧と中心角の関係のグラフを作成します。

# ラベル用の文字列を作成
def_label <- paste0(
  "list(", 
  "phi == frac(1 + sqrt(5), 2), ", 
  "a + b == 2 * pi * r, ", 
  "a + b : a == 1 : frac(1, phi), ", 
  "a + b : b == 1 : 1 - frac(1, phi), ", 
  "alpha == frac(a, r), ", 
  "beta == frac(b, r), ", 
  "theta*degree == frac(180*degree, pi) ~ theta", 
  ")"
)
var_label <- paste0(
  "list(", 
  "r == ", round(r, digits = 2), ", ", 
  "a == ", round(a, digits = 2), ", ", 
  "b == ", round(b, digits = 2), ", ", 
  "alpha == ", round(alpha/pi, digits = 2), " * pi, ", 
  "beta == ", round(beta/pi, digits = 2), " * pi, ", 
  "alpha*degree == ", round(alpha/pi*180, digits = 1), "*degree, ", 
  "beta*degree == ", round(beta/pi*180, digits = 1), "*degree", 
  ")"
)

# ラベル位置の調整値を指定
d_angle <- 0.35
d_arc   <- 0.15

# 円周上で黄金角を作図
add_size <- 0.5
ggplot() + 
  geom_segment(data = angle_axis_df, 
               mapping = aes(x = 0, y = 0, xend = x, yend = y, group = factor(i)), 
               linetype = "dotted") + # 角度目盛線
  geom_text(data = angle_axis_df, 
            mapping = aes(x = x, y = y, label = rad_label, hjust = h, vjust = v), 
            parse = TRUE, size = 4) + # 角度目盛ラベル
  geom_segment(data = radius_df, 
               mapping = aes(x = 0, y = 0, xend = x, yend = y), 
               linewidth = 1) + # 半径線
  geom_path(data = angle_arc_df, 
            mapping = aes(x = x, y = y, group = arc_label), 
            linewidth = 0.5) + # 角マーク
  geom_path(data = circle_df, 
            mapping = aes(x = x, y = y, color = arc_label), 
            linewidth = 1) + # 円周
  geom_text(data = radius_df, 
            mapping = aes(x = 0.5*x, y = 0.5*y, label = edge_label, vjust = -0.5), 
            parse = TRUE, size = 6) + # 半径ラベル
  geom_text(data = angle_label_df, 
            mapping = aes(x = d_angle*x, y = d_angle*y, label = angle_label), 
            parse = TRUE, size = 6) + # 角ラベル
  geom_text(data = angle_label_df, 
            mapping = aes(x = (r-d_arc)*x, y = (r-d_arc)*y, label = arc_label), 
            parse = TRUE, size = 6) + # 弧ラベル
  geom_label(mapping = aes(x = -Inf, y = Inf, label = var_label), 
             parse = TRUE, hjust = 0, vjust = 1, label.r = unit(0, units = "pt"), 
             alpha = 0.5) + # 変数ラベル
  coord_fixed(ratio = 1, 
              xlim = c(-r-add_size, r+add_size), 
              ylim = c(-r-add_size, r+add_size)) + 
  labs(title = "golden angle", 
       subtitle = parse(text = def_label), 
       color = "arc", 
       x = expression(x == r ~ cos~theta), 
       y = expression(y == r ~ sin~theta))

円周の黄金比と黄金角の関係

 円周を黄金角の比  \varphi : 1 = 1 : \frac{1}{\varphi} で分割した円弧を  a a を除いた円弧を  b、また弧  a, b それぞれの中心角を  \alpha, \beta とします。このとき、小さい方の角度  \beta を黄金角と呼びます。
 円弧と中心角(ラジアン) は  \alpha = \frac{a}{r}, \beta = \frac{b}{r} の関係です。

 円周と円弧の関係を線分で表します。

・作図コード(クリックで展開)

 円周ラベルと弧ラベルの描画用のデータフレームを作成します。

# 辺ラベルの座標を作成
d <- 0.1
edge_label_df <- tibble::tibble(
  x = c(0.5*ab, 0.5*a, a+0.5*b), 
  h = 0.5, 
  v = c(-0.5, 1.5, 1.5), 
  edge_label = c("ab", "a", "b")
)
edge_label_df
# A tibble: 3 × 4
      x     h     v edge_label
  <dbl> <dbl> <dbl> <chr>     
1  6.28   0.5  -0.5 ab        
2  3.88   0.5   1.5 a         
3 10.2    0.5   1.5 b    

 円周と円弧を直線にした各線分の中点にラベルを表示することにします。

 角度目盛に対応する目盛の描画用のデータフレームを作成します。

# 半円における目盛数(分母の値)を指定
denom <- 6

# 角度目盛線の座標を作成
angle_axis_df <- tibble::tibble(
  i     = 0:(2*denom),    # 目盛位置番号(分子の値)
  theta = i / denom * pi, # 目盛値
  x     = r * theta, 
  rad_label = paste0("frac(", i, ", ", denom, ") ~ pi * r") # 角度ラベル
)
angle_axis_df
# A tibble: 13 × 4
       i theta     x rad_label           
   <int> <dbl> <dbl> <chr>               
 1     0 0      0    frac(0, 6) ~ pi * r 
 2     1 0.524  1.05 frac(1, 6) ~ pi * r 
 3     2 1.05   2.09 frac(2, 6) ~ pi * r 
 4     3 1.57   3.14 frac(3, 6) ~ pi * r 
 5     4 2.09   4.19 frac(4, 6) ~ pi * r 
 6     5 2.62   5.24 frac(5, 6) ~ pi * r 
 7     6 3.14   6.28 frac(6, 6) ~ pi * r 
 8     7 3.67   7.33 frac(7, 6) ~ pi * r 
 9     8 4.19   8.38 frac(8, 6) ~ pi * r 
10     9 4.71   9.42 frac(9, 6) ~ pi * r 
11    10 5.24  10.5  frac(10, 6) ~ pi * r
12    11 5.76  11.5  frac(11, 6) ~ pi * r
13    12 6.28  12.6  frac(12, 6) ~ pi * r

 角度目盛のときと同様に座標を計算します。ただし、 \theta = 2 \pi を含み、x軸の値は  x = r \theta です。

 円周と円弧の関係のグラフを作成します。

# ラベル用の文字列を作成
def_label <- paste0(
  "list(", 
  "phi == frac(1 + sqrt(5), 2), ", 
  "a + b == 2 * pi * r, ", 
  "a + b : a == phi : 1, ", 
  "a + b : b == phi : phi - 1, ", 
  ")"
)
var_label <- paste0(
  "list(", 
  "r == ", round(r, digits = 2), ", ", 
  "a + b == ", round(ab, digits = 2), ", ", 
  "a == ", round(a, digits = 2), ", ", 
  "b == ", round(b, digits = 2), 
  ")"
)

# 線分上で黄金比を作図
ggplot() + 
  geom_vline(data = angle_axis_df, 
             mapping = aes(xintercept = x), 
             linetype = "dotted") + # 角度目盛線
  geom_line(data = circle_df, 
            mapping = aes(x = r*theta, y = 0, color = arc_label), 
            linewidth = 1) + # 円周
  geom_text(data = edge_label_df, 
            mapping = aes(x = x, y = 0, label = edge_label, hjust = h, vjust = v), 
            parse = TRUE, size = 6) + # 辺ラベル
  geom_vline(mapping = aes(xintercept = a), 
             linetype = "dashed") + # 分割位置線
  geom_text(mapping = aes(x = a, y = Inf), 
            label = expression(frac(2*pi*r, phi)), parse = TRUE, vjust = 1) + # 分割位置ラベル
  geom_label(mapping = aes(x = -Inf, y = Inf, label = var_label), 
             parse = TRUE, hjust = 0, vjust = 1, label.r = unit(0, units = "pt"), 
             alpha = 0.5) + # 変数ラベル
  scale_x_continuous(sec.axis = sec_axis(trans = ~., 
                                         breaks = angle_axis_df[["x"]], 
                                         labels = parse(text = angle_axis_df[["rad_label"]]))) + # 角度目盛ラベル
  scale_y_continuous(breaks = NULL) + 
  labs(title = "golden ratio", 
       subtitle = parse(text = def_label), 
       color = "arc", 
       x = expression(x == r * theta), 
         y = "")

円周と黄金比の関係

 円周の長さが  a + b = 2 \pi r なのを確認できます。

 黄金角の定義より、円周と円弧の長さは、次の関係です。

 \displaystyle
\begin{aligned}
(a + b) : a
   &= \varphi : 1
\\
   &= 1 : \frac{1}{\varphi}
\end{aligned}

 また、 b = (a + b) - a なので、次の関係も成り立ちます。

 \displaystyle
\begin{aligned}
(a + b) : b
   &= \varphi : \varphi - 1
\\
   &= 1 : 1 - \frac{1}{\varphi}
\end{aligned}

 よって、円弧の長さは、次の式で求まります。

 \displaystyle
\begin{aligned}
a  &= (a + b) \frac{1}{\varphi}
\\
   &= 2 \pi r \frac{1}{\varphi}
\\
b  &= (a + b) \left(
          1 - \frac{1}{\varphi}
       \right)
\\
   &= 2 \pi r \left(
          1 - \frac{1}{\varphi}
      \right)
\end{aligned}

 同様に、中心角は、次の式で求まります。

 \displaystyle
\begin{aligned}
\alpha
   &= 2 \pi \frac{1}{\varphi}
\\
   &= \frac{a}{r}
\\
\beta
   &= 2 \pi \left(
          1 - \frac{1}{\varphi}
      \right)
\\
   &= \frac{b}{r}
\end{aligned}

  \alpha, \beta はラジアンです。

 円弧と同じ長さの円周を更に黄金比で分割します。

・作図コード(クリックで展開)

 円周と半径を再設定します。

# 円周の長さを設定
ab <- a

# 半径を計算
r <- ab / (2 * pi)
ab; r
[1] 7.766444
[1] 1.236068

 円周 aba に置き換えて、半径  r = \frac{a + b}{2 \pi} を計算します。

 先ほどのコードで作図します。

黄金比と黄金角の関係

 黄金比により分割を繰り返しても黄金比や黄金角の関係が成り立つのを確認できます。また、円のサイズ(円周の長さ)に関わらず成り立つのも確認できます。

黄金角の計算式

 黄金数を用いない黄金角の計算式を導出します。

 黄金数は、次の式で定義されます。

 \displaystyle
\varphi
    = \frac{1 + \sqrt{5}}{2}

 また、黄金数の逆数は、次のように変形できます(有理化します)。

 \displaystyle
\begin{align}
\frac{1}{\varphi}
   &= \frac{2}{1 + \sqrt{5}}
\\
   &= \frac{2}{1 + \sqrt{5}}
      \frac{1 - \sqrt{5}}{1 - \sqrt{5}}
\\
   &= \frac{2 - 2 \sqrt{5}}{1 - 5}
\\
   &= - \frac{1 - \sqrt{5}}{2}
\tag{1}
\end{align}

 黄金数の変形については「黄金数の性質の導出」を参照してください。

 半径が  r の円周  2 \pi r を黄金比  1 : \frac{1}{\varphi} で円弧  a, b に分割したとき、円周  a + b と短い方の円弧  b は次の比になるのでした。

 \displaystyle
\begin{aligned}
(a + b) : b
   &= 1 : 1 - \frac{1}{\varphi}
\end{aligned}

 またこの式から、円弧  b と中心角  \beta は、次の式で表わせるのでした。

 \displaystyle
\begin{aligned}
b
   &= 2 \pi r \left(
          1 - \frac{1}{\varphi}
      \right)
\\
\beta
   &= 2 \pi \left(
          1 - \frac{1}{\varphi}
      \right)
\end{aligned}

 黄金角  \beta に関する比に式(1)を代入して整理します。

 \displaystyle
\begin{aligned}
1 - \frac{1}{\varphi}
   &= \frac{2}{2} + \frac{1 - \sqrt{5}}{2}
\\
   &= \frac{3 - \sqrt{5}}{2}
\end{aligned}

 この式を黄金角の式に代入します。

 \displaystyle
\begin{aligned}
\beta
   &= 2 \pi
      \frac{3 - \sqrt{5}}{2}
\\
   &= (3 - \sqrt{5}) \pi
\end{aligned}

 黄金角の計算式が得られました。 \beta は弧度法の角度(ラジアン)です。

 度数法の角度の場合は、次の式になります。

 \displaystyle
\beta^{\circ}
    = (3 - \sqrt{5}) 180^{\circ}

 度数法と度数法の角度は  \beta = \frac{2 \pi}{360^{\circ}} \beta^{\circ} の関係です。

 この記事では、黄金角の定義を確認しました。次の記事では、性質を確認します。

おわりに

 フェルマーの螺旋の補足説明用に黄金角とは何ぞやを軽く書くだけのつもりだったのですが、黄金角の性質の図を作ってたら点ごとに中心から距離をとってあげればフェルマーの螺旋になるじゃないかと気付きました。上手く橋渡しする内容になったと思うので諸々の記事をあわせて読んでほしいのですが、思いの外文量があったので分割してアップしていくことにしました。
 それもこれも定義がややこしい所為です。図を見ながら「円周を分割した弧の短い方に対応する角度」と指させばそうなんですが、「円周と円弧」「円弧aとb」「円弧と中心角」「ラジアンと角度」「円の状態と線分の状態」などの対応を1つずつ整理しながら進める必要がありました。単位円で考えれば単純化(例えば弧と角度が一致するとか)できたのですが、単純化することで区別がつかなくなるのもややこしさの原因になり得ますしね。
 「❝わかる❞ために❝わける❞のは基本だから」と今日読んだ本でも言ってました。この本の内容についても近いうちに書く予定です。

 2023年10月29日は、つんく♂さんの55歳のお誕生日です。

 私がハロプロにハマったときにはもう直接プロデュースはしてなかったので、特別深い思い入れがあるわけではないのですが、はっきりとした一番古い記憶だとドラマ「ムコ殿」でのプロデューサー役?事務所の社長役?の印象です。なので一番最初に好きになった曲は「ひとりぼっちのハブラシ」だったかと思います。これからも素敵な曲を楽しみにしています。おめでとうございます!

【次の内容】

www.anarchive-beta.com