フォーゲル螺旋のパッキング密度

# はじめに

https://www.kyushu-u.ac.jp/f/59085/24_10_04.pdf

九大のmath for industryの記事で、パッキング密度という言葉を見かけて、気になったので、コーディング。
「ヒマワリの種の配列は、ほぼ等間隔」という話から、黄金角以外はどうなの？という例(ためし)。

# フェルマー螺旋とフォーゲル螺旋の違い

ざっくり同じ式です。
フェルマー螺旋は連続した曲線であり、
フォーゲル螺旋は離散した点群です。

| 特徴 | Fermat Spiral | Vogel Spiral |
|:-:|:-:|:-:|
|数式|r = √θ|r = c√n<br>θ = nφ|
|用途|曲線(連続)|点群(離散)|
|関係性|螺旋|葉、種配置|
|回転角|任意角度|黄金角|

```math
黄金角\phi=\frac{2\pi}{\phi^2}\approx 2.39996 rad
```

# フェルマー螺旋 fermat's spiral

```java
void setup() {
  size(800, 800);
  float r,x,y;  
  translate(width/2, height/2);
  
  beginShape();
  for (float th = 0; th < 100; th += 0.1) {
    r = 50 * sqrt(th);
    x = r * cos(th);
    y = r * sin(th);
    vertex(x, y);
  }
  endShape();
}
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/2453a029-8922-419a-9894-6952308739cb.png)

# フォーゲル螺旋 Vogel Spiral

フェルマー螺旋は連続した曲線で、その螺旋上の飛び飛びの点がフォーゲル螺旋です。

```java
void setup() {
  size(800, 800);
  float r,x,y;
  translate(width/2, height/2);
  fill(0);
  
  for (float th = 0; th < 10000; th += radians(137.50776405)) {
    r = 5 * sqrt(th);
    x = r * cos(th);
    y = r * sin(th);
    circle(x, y,8);
  }
}
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/9a35ee3b-db2d-4501-ab2e-b006e0d7f59e.png)

# フォーゲル螺旋 Vogel Spiralの角度を変更してみる

一般的には黄金角を使うところを、他の数字を使う。九大の資料を参考にした。
画像の左側は、適当にレインボーカラーに着色した点群。
画像の右側は、適当に密度を可視化した。上下が逆。

↓黄金角 : 137.50776405
![GRA.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/7da3e3bf-74bf-425f-9ea5-849362f7d8ea.png)

↓白銀角 : 210.88311757
![SRA.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/c130ed8f-3b1a-4c49-951f-c2ff707707f1.png)

↓青銅角 : 250.00077042
![ARA.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/6cf1f2c0-f9c7-4430-98a0-6431ef42b7fc.png)

↓(9+√221)/10 : 209.18171236
![RRA.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/a8e5b3ce-5f58-4ea1-9e86-86549e15fe3d.png)

↓PI : 245.40844096
![PIA.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/cf2ac80f-c74b-4edb-af8b-7527a9f56da6.png)

↓e : 227.56340117
![E_A.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/bf728973-dcf0-47c4-8bc4-c4e0e22456fd.png)

```java:VogelSpiral.pde
float GRA = 137.50776405; //Golden Ratio     1.6180339887
float SRA = 210.88311755; //Silver Ratio     2.4142135623
float RRA = 209.18171236; //(9+sqrt(221))/10 2.3866068747
float PIA = 245.40844096; //PI 3.1415926536
float E_A = 227.56340117; //e  2.7182818285
float ARA = 250.00077042; //3.30277563773199

int total = 13000;
float angle = radians(GRA);
float c = 5;

int gridSize = 10;  // 密度マップのセルサイズ
int cols, rows;
int[][] density;

void setup() {
  size(1600, 800);
  colorMode(HSB, 255);
  background(0);
  translate(400, 400);
  
  cols = width / gridSize;
  rows = height / gridSize;
  density = new int[cols][rows];

  for (int n = 0; n < total; n++) {
    float a = n * angle;
    float r = c * sqrt(n);
    float x = r * cos(a);
    float y = r * sin(a);

    float hue = map(n, 0, total, 0, 255);

    noStroke();
    fill(hue, 200, 255);
    ellipse(x, y, 4, 4);

    int cx = int((x + 400) / gridSize);
    int cy = int((y + 400) / gridSize);
    
    if (cx >= 0 && cx < cols && cy >= 0 && cy < rows) {
      density[cx][cy]++;
    }
  }
  
  // 密度最大値を調べる
  int maxDensity = 0;
  for (int i = 0; i < cols; i++) {
    for (int j = 0; j < rows; j++) {
      if (density[i][j] > maxDensity) {
        maxDensity = density[i][j];
      }
    }
  }

  // 可視化（ヒートマップ）
  noStroke();
  translate(400, -400);
  for (int i = 0; i < cols; i++) {
    for (int j = 0; j < rows; j++) {
      float hue = map(density[i][j], 0, maxDensity, 0, 255);
      fill(hue, 255, 255);
      rect(i * gridSize, j * gridSize, gridSize, gridSize);
    }
  }
  //save("GRA.png");
}
```

# 黄金角の精度問題（137.5では足りない？）

黄金角は137.50776405...と永遠に綴らなる数字ですが、他の資料では、137.5で計算していることがあります。その次の桁が０であることから、これで十分と思いがちですが、ヒマワリの種の配列で点をたくさん描いていくと、誤差が積もり積もって、螺旋模様が消えてしまいます。
以下は、精度を上げていった場合の、見た目の変化を示しています。
点の数は80000まで増やし、全体も縮小しています。
中央ほど、精度が足りなくても、問題ありませんが、
周辺に行くほど、誤差の影響が大きくなります。


ただし、floatの精度の都合上、
137.50777
(137.50776を代入するとこうなる)
が限界でした。これ以上は、doubleに対応した算術関数を使うなどの対策が必要です。
というわけで、
137.5ではなく、めんどくさがらずに、137.50776を使いましょう。

↓137.5
![GRA137_5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/8350ad09-959f-4d3a-b991-deddb085a97c.png)
↓137.507
![GRA137_507.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/91102369-171d-419e-b781-b43a06f3073e.png)
↓137.5077
![GRA137_5077.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/092ed055-103a-4f36-b13e-93f47e776fc0.png)
↓137.50777
![GRA137_50776.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/486a0a10-f74f-4d67-b04a-7962a8fd5d1c.png)


↓いろいろ比較してるサイト。

https://www.codeproject.com/Articles/1221341/The-Vogel-Spiral-Phenomenon