iGeo Tutorials

チュートリアル

ニュートンの運動法則とパーティクル

位置、速度、加速度、力

古典力学では、物の動きは、位置、速度、加速度、力の各々のベクトル変数と、質量のスカラ変数で記述されます。速度は位置の時間微分で与えられ、加速度は速度の時間微分で与えられ、位置の２階微分ともなります。

一方、力と質量と加速度は、ニュートンの運動方程式によって関係づけられています。

ニュートンの運動法則に従ってパーティクルの動きをシミュレートするエージェント・クラスを実装するには、まず位置、速度、加速度、力、質量の変数をフィールドとして追加します。以下のコードでは質量はdouble型、他のフィールドはIVec型としています。通常、パーティクルのシミュレーションは初期位置と初速度を与えて実行するので、コンストラクタの引数に位置と速度の２つの変数をとります。

import processing.opengl.*;
import igeo.*;

class MyParticle extends IAgent{
  IVec pos;
  IVec vel;
  IVec acc;
  IVec frc;
  double mass = 1.0; // default mass
  
  MyParticle(IVec p, IVec v){
    pos = p;
    vel = v;
  }
}

位置と速度はパーティクルの内部状態と言えますが、力は外部から与えられるものです。シミュレーションのタイム・フレームごとに、パーティクルは力を受け、それに応じて速度と位置が更新されます。一つのパーティクルが同時に複数の力受けるとき、それは、それらのベクトル和の力を一回受けることと等しくみなせます。

以下の図にベクトル和の例を示します。

パーティクル・エージェントが力を外部または内部から受けるための push()メソッドを以下のように実装します。

import processing.opengl.*;
import igeo.*;

class MyParticle extends IAgent{
  IVec pos, vel, acc, frc;
  double mass = 1.0; // default mass
  
  MyParticle(IVec p, IVec v){
    pos = p;
    vel = v;
    frc = new IVec();
  }

  void push(IVec f){
    frc.add(f);
  }
}

複数の力の各々の力ベクトルでなく、合算された力のベクトル和が一つあれば加速度、速度、位置は計算できます。加速度、速度、位置は数学的には以下のように計算されます。

パーティクル・クラスの実装

上に示された式をエージェントの中で実装するために、 update()メソッドを以下のように書きます。

import processing.opengl.*;
import igeo.*;

class MyParticle extends IAgent{
  final double dt = IConfig.updateRate; //second per frame
  IVec pos, vel, acc, frc;
  double mass = 1.0; // default mass
  
  MyParticle(IVec p, IVec v){
    pos = p;
    vel = v;
    frc = new IVec();
  }

  void push(IVec f){
    frc.add(f);
  }

  void update(){
    acc = frc.div(mass);
    vel.add(acc.mul(dt));
    pos.add(vel.cp().mul(dt)); //vel itself can't be changed
    frc.zero(); //reset frc by setting zero
  }
}

変数dtはエージェントが更新されるタイム・フレーム間の秒数です。(通常は0.03-0.1秒程度) この変数は定数なのでfinalというキーワードが付加されます。 iGeoのデフォルトの更新秒数はIConfig.updateRateで設定されています。 iGeoシステムの更新秒数はIG.updateRate(double)メソッドやその別名短縮メソッドIG.rate(double)やIG.speed(double)で設定することができます。アップデート・メソッドでは、加速度accが力frcをmassで割ることによって計算されています。それから速度velは、微小更新時間分の加速度の累積を加算することで更新されます。具体的にはaccにdtを掛けたものを足します。同様に位置を更新するには、速度velに微小時間dtかけたものを足しています。ただし、速度は元の値を保持する必要があるため、cp()メソッドで複製されてから掛け算されています。そしてアップデート・メソッドの最後では力frcをゼロにします。

パーティクル・エージェントを用いるとき、だいたいは何らかの力をパーティクルに作用し、なおかつ、その中でなんらかの幾何学オブジェクトを生成します。以下の例では、パーティクル・エージェント自身が自分にY軸の下の方向に力を与え、またそれぞれの時点での位置に点を生成して行きます。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  IG.duration(200);
  new MyParticle(IG.v(0,0,0), IG.v(10,20,0));
}

class MyParticle extends IAgent{
  final double dt = IConfig.updateRate;//second per frame
  IVec pos, vel, acc, frc;
  double mass = 1.0; // default mass
  
  MyParticle(IVec p, IVec v){
    pos = p;
    vel = v;
    frc = new IVec();
  }

  void push(IVec f){
    frc.add(f);
  }

  void update(){
    push(IG.v(0,-10,0));
    new IPoint(pos.cp());
    
    acc = frc.div(mass);
    vel.add(acc.mul(dt));
    pos.add(vel.cp().mul(dt)); //vel itself can't be changed
    frc.zero(); //reset frc by setting zero
  }
}

以下の例ではコードの構造化のために、パーティクルの動きをシミュレートする部分と、幾何学オブジェクトを生成する部分を分離し、前者をMyParticleBaseに入れ、後者をMyParticleBaseを継承するMyParticleに入れました。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  IG.duration(200);
  new MyParticle(IG.v(0,0,0), IG.v(10,20,0));
}

class MyParticle extends MyParticleBase{
  MyParticle(IVec p, IVec v){ super(p,v); }

  void update(){
    push(IG.v(0,-10,0));
    new IPoint(pos.cp());
    super.update();
  }
}

class MyParticleBase extends IAgent{
  final double dt = IConfig.updateRate;//second per frame
  IVec pos, vel, acc, frc;
  double mass = 1.0; // default mass
  
  MyParticleBase(IVec p, IVec v){
    pos = p;
    vel = v;
    frc = new IVec();
  }

  void push(IVec f){ frc.add(f); }

  void update(){
    acc = frc.div(mass);
    vel.add(acc.mul(dt));
    pos.add(vel.cp().mul(dt)); //vel itself can't be changed
    frc.zero(); //reset frc by setting zero
  }
}

IParticleクラスの利用

iGeoライブラリには、パーティクルの動作をシミュレートする IParticleクラスが用意されています。このクラスは上記の例で記したような実装されています。一つ上のコードは、 MyParticleBaseをIParticleで置き換えて、以下のように書き直すことができます。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  IG.duration(200);
  new MyParticle(IG.v(0,0,0), IG.v(10,20,0));
}

class MyParticle extends IParticle{
  MyParticle(IVec p, IVec v){ super(p,v); }

  void update(){
    push(IG.v(0,-10,0));
    new IPoint(pos().cp());
  }
}

また、IParticleクラスを継承して独自の子クラスを作成するときには、アップデート・メソッドに super.update(); と記述する必要はありません。

以下はIParticleを継承して独自の振る舞いをするパーティクル・エージェントのもう一つの例です。アップデート・メソッドでランダムな力を自身に与えています。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  new MyParticle(IG.v(0,0,0), IG.v(0,0,0));
  new MyParticle(IG.v(0,0,0), IG.v(0,0,0));
  new MyParticle(IG.v(0,0,0), IG.v(0,0,0));
}

class MyParticle extends IParticle{
  MyParticle(IVec p, IVec v){ super(p,v); }

  void update(){
    IVec force = IRand.pt(-100,-100,0,100,100,0);
    push(force);
    new IPoint(pos().cp());
  }
}

もう一つの例は、パーティクルの速度に応じて計算した力ベクトルを自身に適用しています。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  new MyParticle(IG.v(0,0,0), IG.v(10,0,0));
}

class MyParticle extends IParticle{
  MyParticle(IVec p, IVec v){ super(p,v); }

  void update(){
    IVec force = vel().cp().mul(2).rot(PI/2);
    push(force);
    new IPoint(pos().cp());
  }
}

パーティクル間のインタラクション

以下のコードはパーティクル同士がインタラクト・メソッドを通じてインタラクトする例です。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  for(int i=0; i < 10; i++){
    new MyParticle(IG.v(0,i*10,0), IG.v(10,0,0)).clr(0.1*i,0,0);
  }
}

class MyParticle extends IParticle{
  MyParticle(IVec p, IVec v){ super(p,v); }

  void interact(ArrayList < IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof MyParticle){
        MyParticle p = (MyParticle)agents.get(i);
        if(p != this){
          if(p.dist(this) < 20){ //distance threshold
            IVec force = p.pos().dif(pos());
            force.mul(0.01);
            push(force);
          }
        }
      }
    }
  }
  
  void update(){
    new IPoint(pos().cp()).clr(clr());
  }
}

インタラクト・メソッドでパーティクルは、他のパーティクルの位置をチェックし、もしそれがある距離より近ければ、相手のパーティクルからの差ベクトルを計算し、それにある値をかけ、自身に力ベクトルとして適用しています。その結果、パーティクルはお互いが引き合っているような振る舞いをみせます。

以下の例は、上記にある例の、速度に応じた力で回転するパーティクル・エージェントに、ひとつ前の例と同様の、引き合う動作をインタラクト・メソッドに定義し、振る舞いの変化を見るためのものです。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  for(int i=0; i < 10; i++){
    new MyParticle(IRand.pt(0,0,100,100), IG.v(10,0,0)).clr(IRand.clr());
  }
}

class MyParticle extends IParticle{
  MyParticle(IVec p, IVec v){ super(p,v); }

  void interact(ArrayList < IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof MyParticle){
        MyParticle p = (MyParticle)agents.get(i);
        if(p != this){
          if(p.dist(this) < 20){
            IVec force = p.pos().dif(pos());
            force.mul(1.2);
            push(force);
          }
        }
      }
    }
  }
  
  void update(){
    IVec force = vel().cp().mul(2).rot(PI/2);
    push(force);
    new IPoint(pos().cp()).clr(clr());
  }
}

次のコードはランダム・ベクトルをランダム・ウォークの動作と引き合うインタラクト・メソッドを組み泡あせたものです。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  for(int i=0; i < 30; i++){
    new MyParticle(IRand.pt(0,0,100,100), IG.v(10,0,0)).clr(IRand.clr());
  }
}

class MyParticle extends IParticle{
  MyParticle(IVec p, IVec v){ super(p,v); }

  void interact(ArrayList < IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof MyParticle){
        MyParticle p = (MyParticle)agents.get(i);
        if(p != this){
          if(p.dist(this) < 20){
            IVec force = p.pos().dif(pos());
            force.mul(2.0);
            push(force);
          }
        }
      }
    }
  }
  
  void update(){
    IVec force = IRand.pt(-100,-100,0,100,100,0);
    push(force);
    new IPoint(pos().cp()).clr(clr());
  }
}

パーティクルに力を及ぼす力の場エージェント

これまでの例では、パーティクルに及ぼされる力はパーティクル自身が与えていました。しかし、重力のように外部から与えられる力も多くあります。以下ではこのような、パーティクルに力を与える力の場エージェントを実装します。力の場エージェントはインタラクト・メソッドのforループでパーティクルをみつけだして、力を適用します。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  IG.duration(200);
  new MyParticle(IG.v(0,0,0), IG.v(10,20,0));
  new MyGravity();
}

class MyParticle extends IParticle{
  MyParticle(IVec p, IVec v){ super(p,v); }
  void update(){ 
    new IPoint(pos().cp()); 
  }
}

class MyGravity extends IAgent{

  void interact(IDynamics agent){
    if(agent instanceof MyParticle){
      MyParticle particle = (MyParticle)agent;
      particle.push(IG.v(0,-10,0));
    }
  }
}

上記の例では、短縮版のインタラクトメソッド
void interact(IDynamics agent)
を用いています。インタラクト・メソッドの中のこのif条件分岐
if( agent instanceof MyParticle )
では、該当エージェントがMyParticleのインスタンスであるかどうかをチェックしています。もしそうであれば、IDynamics型の変数agentは MyParticle型の変数particleにキャストにより変換されたのち、 push()メソッドによって力が適用されます。 interact()メソッドの詳細については該当チュートリアル参照してください。

次のコードは上のコードを少し変更して、 MyGravityクラスにIVec型のフィールドgravityを加えて、コンストラクタで重力の方向と大きさを設定可能にしてあります。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  IG.duration(200);
  new MyParticle(IG.v(0,0,0), IG.v(10,20,0));
  new MyGravity(IG.v(0,-10,0));
}

class MyParticle extends IParticle{
  MyParticle(IVec p, IVec v){ super(p,v); }
  void update(){ new IPoint(pos().cp()); }
}

class MyGravity extends IAgent{
  IVec gravity;
  
  MyGravity(IVec g){ gravity=g; }

  void interact(IDynamics agent){
    if(agent instanceof MyParticle){
      MyParticle particle = (MyParticle)agent;
      particle.push(gravity);
    }
  }
}

以下のコードでは、初期位置と初速度の異なるパーティクルに重力エージェントを適応しています。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  IG.duration(200);
  int num = 40;
  for(int i=0; i < num; i++){
    double a = i*2*PI/num;
    new MyParticle(IG.v(20,0,0).rot(a),IG.v(10,0,i).rot(a));
  }
  new MyGravity(IG.v(0,0,-10));
}

class MyParticle extends IParticle{
  MyParticle(IVec p, IVec v){ super(p,v); }
  void update(){ new IPoint(pos().cp()); }
}

class MyGravity extends IAgent{
  IVec gravity;
  
  MyGravity(IVec g){ gravity=g; }

  void interact(IDynamics agent){
    if(agent instanceof MyParticle){
      MyParticle particle = (MyParticle)agent;
      particle.push(gravity);
    }
  }
}

パーティクルを跳ね返す平面エージェント

以下では、XY平面上でパーティクルを跳ね返させるエージェントを作成します。 BouncePlaneと名付けられたクラスには、フィールドzが与えられ、平面のz軸上での位置を決めます。インタラクト・メソッドでは、パーティクルのz座標ががフィールドzよりも低い位置にいるかどうかチェックされ、もしそうであればまずその位置まで押し上げられ、それからパーティクルの速度ベクトルを z軸方向へ反転させます。あるベクトル方向への反転はrefメソッドを用います。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  IG.duration(200);
  int num = 40;
  for(int i=0; i < num; i++){
    double a = i*2*PI/num;
    new MyParticle(IG.v(20,0,0).rot(a), IG.v(10,0,i).rot(a));
  }
  new MyGravity(IG.v(0,0,-10));
  new BouncePlane(-5);
}

class MyParticle extends IParticle{
  MyParticle(IVec p, IVec v){ super(p,v); }
  void update(){ new IPoint(pos().cp()); }
}

class MyGravity extends IAgent{
  IVec gravity;
  MyGravity(IVec g){ gravity=g; }
  void interact(IDynamics agent){
    if(agent instanceof MyParticle){
      MyParticle particle = (MyParticle)agent;
      particle.push(gravity);
    }
  }
}

class BouncePlane extends IAgent{
  double z;
  
  BouncePlane(double zval){ z=zval; }

  void interact(IDynamics agent){
    if(agent instanceof MyParticle){
      MyParticle particle = (MyParticle)agent;
      if(particle.pos().z() < z){
        particle.pos().add(0,0,z - particle.pos().z());
        particle.vel().ref(IG.zaxis);
      }
    }
  }
}

パーティクルの軌跡としての幾何学オブジェクト生成

パーティクルの軌跡上になんらかの幾何学オブジェクトを作成する場合、アップデート・メソッドの中で幾何学オブジェクトを生成します。もし１フレーム前の位置から現在の位置に線を引く場合は、前の位置をprevPosのような変数に覚えておかなければなりません。また、一番最初のフレームでは以前の位置が存在しないため、 prevPosにはnull値になっており、そのケースを除外するためにif条件分岐がなされています。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  IG.duration(200);
  int num = 40;
  for(int i=0; i < num; i++){
    double a = i*2*PI/num;
    new MyParticle(IG.v(20,0,0).rot(a), IG.v(10,0,i).rot(a));
  }
  new MyGravity(IG.v(0,0,-10));
  new BouncePlane(-5);
}

class MyParticle extends IParticle{
  IVec prevPos;
  MyParticle(IVec p, IVec v){ super(p,v); }
  void update(){
    IVec curPos = pos().cp();
    if(prevPos != null){
      new ICurve(prevPos, curPos).clr(0);
    }
    prevPos = curPos;
  }
}

class MyGravity extends IAgent{
  IVec gravity;
  MyGravity(IVec g){ gravity=g; }
  void interact(IDynamics agent){
    if(agent instanceof MyParticle){
      MyParticle particle = (MyParticle)agent;
      particle.push(gravity);
    }
  }
}

class BouncePlane extends IAgent{
  double z;
  
  BouncePlane(double zval){ z=zval; }

  void interact(IDynamics agent){
    if(agent instanceof MyParticle){
      MyParticle particle = (MyParticle)agent;
      if(particle.pos().z() < z){
        particle.pos().add(0,0,z - particle.pos().z());
        particle.vel().ref(IG.zaxis);
      }
    }
  }
}

上の例で先端に残っている点は、hide()メソッドで隠すことができます。

軌跡上に幾何学オブジェクトを生成するもう一つの方法は、幾何学オブジェクトを生成ためだけの別のエージェント・クラスを作成することです。以下の例では GeometrAgentクラスが追加され、その内部に２つのパーティクルとそれぞれの以前の位置をフィールドとして持ちます。コンストラクタでパーティクルのペアが与えられ、アップデート・メソッドにおいて２つのパーティクルとそれぞれの一つ以前の位置の合計４点から矩形の面を生成します。

import processing.opengl.*;
import igeo.*;

void setup(){
  size(480, 360, IG.GL);
  IG.duration(200);
  int num = 40;
  MyParticle[] particle = new MyParticle[num];
  for(int i=0; i < num; i++){
    double a = i*2*PI/num;
    particle[i] =
    new MyParticle(IG.v(20,0,0).rot(a), IG.v(10,0,i).rot(a));
  }
  for(int i=1; i < num; i++){
    new GeometryAgent(particle[i-1], particle[i]);
  }
  new MyGravity(IG.v(0,0,-10));
  new BouncePlane(-5);
}

class MyParticle extends IParticle{
  IVec prevPos;

  MyParticle(IVec p, IVec v){ 
    super(p,v);
    hide(); // hiding the point of IParticle
  }
  void update(){
    IVec curPos = pos().cp();
    if(prevPos != null){
      new ICurve(prevPos, curPos);
    }
    prevPos = curPos;
  }
}

class MyGravity extends IAgent{
  IVec gravity;
  MyGravity(IVec g){ gravity=g; }
  void interact(IDynamics agent){
    if(agent instanceof MyParticle){
      MyParticle particle = (MyParticle)agent;
      particle.push(gravity);
    }
  }
}

class BouncePlane extends IAgent{
  double z;
  
  BouncePlane(double zval){ z=zval; }

  void interact(IDynamics agent){
    if(agent instanceof MyParticle){
      MyParticle particle = (MyParticle)agent;
      if(particle.pos().z() < z){
        particle.pos().add(0,0,z - particle.pos().z());
        particle.vel().ref(IG.zaxis);
      }
    }
  }
}

class GeometryAgent extends IAgent{
  MyParticle particle1, particle2;
  IVec prevPos1, prevPos2;

  GeometryAgent(MyParticle p1, MyParticle p2){
    particle1 = p1;
    particle2 = p2;
  }

  void update(){
    IVec curPos1 = particle1.pos().cp();
    IVec curPos2 = particle2.pos().cp();
    if(prevPos1 != null && prevPos2 != null){
      new ISurface(prevPos1, curPos1, curPos2, prevPos2).clr(clr());
    }
    prevPos1 = curPos1;
    prevPos2 = curPos2;
    // update color
    double r = red() + IRand.get(-0.01,0.01);
    double g = green() + IRand.get(-0.02,0.02);
    clr(r, g, 1.0);
  }
}

(トピック一覧へ戻る)

HOME

FOR PROCESSING

DOWNLOAD

DOCUMENTS

TUTORIALS (Java / Python)

GALLERY

SOURCE CODE(GitHub)

ABOUT/CONTACT