Processingでシミュレーション～多粒子の運動(Mass-Spring-Damper)

# はじめに
Processingを使って多粒子の運動をシミュレーションします。

粒子同士の反発はMass-Spring-Damperを使ってモデリングします。
外力は重力とします。

粒子は質点と考えます。
粒子同士および壁(ウィンドウ境界)との影響が及ぶ半径を影響半径$r_e$とします。
$r_e$範囲内では、外力$\vec{F_{out}}$とSpring-Damper$\vec{F_{sd}}$による力が働きます。
$r_e$範囲外では、外力のみが働きます。

# 外力
外力は重力のみとします。$m$は粒子の質量です。$g$は重力加速度です。

```math
\vec{F_{out}} = (0, mg) \\
```

# Spring-Damperによる力
2粒子の間の力を考えます。それぞれに番号$i,j$をつけます。
Spring-Damper$\vec{F_{sd}}$による力は次の式で与えられます。

```math
\vec{F_{sd}} = \{kd + c v\} \vec{n} \\
d = r_{ei} + r_{ej} - |\vec{r}_j - \vec{r}_i| \\
v = \{\vec{v}_j - \vec{v}_i \} \cdot \vec{n} \\
```
$k$はばね定数、$c$は減衰係数です。$\vec{n}$は影響半径から遠ざかる方向の単位ベクトルです。
$d$は影響半径に侵入した距離です。$v$は2粒子の相対速度の$\vec{n}$に沿った成分です。

# Mass-Spring-Damper
運動方程式を示します。
$r_e$範囲内では次の式になります。

```math
m \ddot{\vec{r}} = \vec{F_{out}} + \vec{F_{sd}} = \\
(0, mg) + \{kd + c v\} \vec{n} \\
```

$r_e$範囲外では次の式になります。

```math
m \ddot{\vec{r}} = \vec{F_{out}}  = (0, mg) \\
```

## 粒子と壁の接近
粒子が壁に接近した状況を考えます。壁が粒子の$r_e$範囲内に入ると
壁から遠ざかる単位ベクトル$\vec{n}$の方向に$\vec{F}_{sd}$が働きます。

![image](https://qiita-image-store.s3.amazonaws.com/0/84205/1e829525-7154-ae00-d6d2-80b353a4faf1.png)

## 粒子同士の接近
2粒子が接近した状況を考えます。粒子に番号$i,j$をつけます。

![image](https://qiita-image-store.s3.amazonaws.com/0/84205/e674bee3-cbe0-e82a-7bf8-65c9403a3a38.png)

粒子の位置ベクトルを$\vec{r}_i, \vec{r}_j$とします。
2粒子の距離を

```math
r = |\vec{r}_j - \vec{r}_i|
```
とします。

$r$が2粒子の$r_e$の和よりも小さいときに接近したと判定します。

```math
r < r_{ej} + r_{ei} \\
0 <(r_{ej} + r_{ei}) -r = r_e
```

$\vec{F_{sd}}$は$r_e$および相対速度$v$に応じて決まります。

$F_{sd}$の方向
粒子$i$については$j$から$i$に向かう方向の単位ベクトル$\vec{n}_{ji}$となり、
粒子$j$についてはその逆方向になります。

```math
\vec{n}_{ji} = \frac{|\vec{r}_i - \vec{r}_j|}{r} \\
\vec{n}_{ij} = -\vec{n}_{ji} \\
```

## 差分化
粒子の位置は次の順番で計算します。

+ 時刻$t$の$\vec{F}(t)$を求める。
+ $\vec{F}(t)$から時刻$t+Δt$の速度$\vec{v}(t+Δt)$を求める。
+ 速度$\vec{v}(t+Δt)$から時刻$t+2Δt$の位置$\vec{r}(t+2Δt)$を求める。

力
$r_e$範囲内

```math
\vec{F}(t) = (0, mg) \\
```

$r_e$範囲外

```math
\vec{F}(t) = (0, mg) + \{kd + c v\} \vec{n} \\
d = r_{ei} + r_{ej} - |\vec{r}_j(t+Δt) - \vec{r}_i(t+Δt)| \\
v = \{\vec{v}_j(t) - \vec{v}_i(t) \} \cdot \vec{n} \\
\vec{n} = \frac{\vec{r}_i(t+Δt) - \vec{r}_j(t+Δt)}{|\vec{r}_i(t+Δt) - \vec{r}_j(t+Δt)|} \\
```

速度

```math
\vec{v}(t+Δt) = \vec{v}(t) + \frac{Δt}{m} \vec{F}(t) \\
```
位置

```math
\vec{r}(t+2Δt) = \vec{r}(t+Δt) + \vec{v}(t+Δt)Δt \\
```

# Processingで可視化
粒子数N=2000の場合のシミュレーションを示します。
各粒子の速度、位置はランダムにしています。

[![IMAGE ALT TEXT HERE](http://img.youtube.com/vi/JGJ-r5fXeKA/0.jpg)](http://www.youtube.com/watch?v=JGJ-r5fXeKA)

```java
ParticleSet set; //<>//
float k = 500.00; // spring constant 
float c = -1.0; // damping coefficient
float g = 9.80; // unit:m/s/s, gravitational acceleration
float fps = 60;
float dt = 0.01; // unit:s, time step
float mpp = 0.01; // unit:m/pixel
float remin = 0.05, remax = remin; // unit:m, radius
int numOfParticles = 2000;
float mmin = 0.1, mmax = mmin; // unit:kg, range of mass
PVector r0min, r0max, v0min, v0max; // unit:m, m/s, range of initial value
float hsv_max = 100;
float color_s = hsv_max, color_v = hsv_max;
int speed = 1;
boolean collision_enb = true;

void setup() {
  size(640, 480);
  frameRate(fps);

  // initial value
  r0max = new PVector(mpp*width-remax, mpp*height-remax);
  r0min = new PVector(remax, remax);
  v0max = new PVector( 1.0, 1.0);
  v0min = new PVector(-1.0,-1.0);

  set = new ParticleSet(numOfParticles);
  colorMode(HSB, hsv_max, hsv_max, hsv_max);   
}

void draw() {
  background(100);
  for (int i=0; i<speed; i++)
    set.step();
  set.draw();
  // Saves each frame as line-000001.png, line-000002.png, etc.
  //saveFrame("line-######.png");  
}

class ParticleSet {
  Particle[] set;

  ParticleSet (int num) {
    set = new Particle[num];
    for (int i=0; i<set.length; i++) {
      float re = random(remin, remax);
      float x = random(r0min.x, r0max.x);
      float y = random(r0min.y, r0max.y);
      float vx = random(v0min.x, v0max.x);
      float vy = random(v0min.y, v0max.y);
      float m = random(mmin, mmax);
      set[i] = new Particle(new PVector(x, y), new PVector(vx, vy), m, re);
    }
  }

  void draw() {
    for (int i=0; i<set.length; i++) {
      set[i].draw();
    }
  }

  void step() {
    calc_f();
    for (int i=0; i<set.length; i++) {
      set[i].step();
    }
  }

  void calc_f() {
    for (int i=0; i<set.length; i++) {
      set[i].Fsd = new PVector();
    }
    for (int i=0; i<set.length; i++) {
      for (int j=0; j<set.length && j < i; j++) {
        PVector rij = PVector.sub(set[j].r, set[i].r);
        float dij = rij.mag();
        float reij = set[i].re + set[j].re;
        float d = reij - dij;

        if (0 < d) {
          PVector nij = rij.normalize().copy();
          PVector Fs = PVector.mult(nij, k*d);
          PVector vij = PVector.sub(set[j].v, set[i].v);
          float v = PVector.dot(nij, vij);
          PVector Fd = PVector.mult(nij, v*c);
          PVector Fsd = PVector.add(Fs, Fd);
          set[i].Fsd = PVector.sub(set[i].Fsd, Fsd);
          set[j].Fsd = PVector.add(set[j].Fsd, Fsd);
        }
      }
    }
  }
}

class Particle {
  PVector F, Fsd;
  PVector v;
  PVector r;
  float m;
  float re;
  float color_h;

  Particle(PVector r0, PVector v0, float m, float re) {
    r = r0;
    v = v0;
    this.m = m;
    this.re = re;
    color_h = random(0, hsv_max);
  }

  void draw() {
    float vmag = abs(v.x)+abs(v.y); color_h = map(vmag, 0, 10, 70, 0);
    stroke(color_h, color_s, color_v); /* HSV */
    fill(color_h, color_s, color_v, 50); /* HSVA */
    ellipse(r.x/mpp, r.y/mpp, re*2/mpp, re*2/mpp);
  }

  void step() {
    F = new PVector(0, m*g); // Fout
    calc_f();
    if (collision_enb)
      F = PVector.add(F, Fsd);
    v = PVector.add(v, PVector.mult(F, dt/m));
    r = PVector.add(r, PVector.mult(v, dt));
  }

  void calc_f() {
    calc_f_wall();
  }

  void calc_f_wall() {
    float d, n;
    // x direction
    if (r.x < re) {
      n =  1.0;
      d = re-r.x;
    } else if (mpp*width < r.x+re) {
      n = -1.0;
      d = r.x+re-mpp*width;
    } else {
      n = 0;
      d = 0;
    }
    if (0 < d) {
      F.x += k*d*n;
      F.x += c*v.x;
    }

    // y direction
    if (r.y < re) {
      n =  1.0;
      d = re-r.y;
    } else if (mpp*height < r.y+re) {
      n = -1.0;
      d = r.y+re-mpp*height;
    } else {
      n = 0;
      d = 0;
    }
    if (0 < d) {
      F.y += k*d*n;
      F.y += c*v.y;
    }
  }
}
```