今回作ってもらう物理シミュレーションの「土台」としてPhystemというもの(オブジェクト)が用意されている。
このファイルではそのPhystemの使い方と、動作原理を述べておこう。
コンピュータの上に物体の運動を再現したい。その為に画面に物体を「描く」のだが、ここでは2次元で直交座標で扱うことにする。
Phystemはその中にjavascriptのcanvasを持っている。canvasの詳しい使い方を知りたい人はjavascriptに関する文献(たとえば、html5.jpのCanvasリファレンスなど、ネット上のリファレンスを読むとよいが、ここでは最低限必要なことだけを記しておく。
コンピュータ上でよく使われる座標系は「ピクセル」(液晶やCRTの輝点一つを1とする単位)になっていて、かつ右がx軸の正、上がy軸の正と決めてある場合が多い(そして原点はたいてい、画面の左上隅である)。
図で示すと↑のような感じになっている。しかしこれは物理に使うには使いにくいので、まず座標系の向きは右がx軸の正、上がy軸 の正とする。原点は任意の場所に置けるようにし、かつ単位系もピクセルとは関係ない自由なスケールにしておく(縦横比はピクセルの比のまま1:1である)。
Phystemの中では座標系のいろいろな値が上のように決められている。たとえばPhystemを表す変数がpsであったなら、ps.wが横幅を、 ps.leftxが左隅x座標、ps.rightxが右隅x座標である(当然であるが、ps.w = ps.rightx - ps.leftxである)。図では原点(0,0)が中にあるが、別に中になくてもよい。ps.canvasがキャンバス、ps.ctxが描画コンテキストなので、これを 使って絵を描くことができる。
プログラムの最初で、
のように書いておくと、横幅20、縦幅14の描画領域がidが「canvas1」←これは適当に変えてよいのcanvasとして作られる。
canvasの場所は、ファイルの中の
と書いた部分である。id=の後の部分が一致していることに注意。
最後の引数の0.8は、このcanvasがブラウザの画面上の横幅の何倍になるかを指定する。
この作り方だと、中央に原点が来る。つまり、x座標は-10から10に、y座標が-7から7になる。これを変更したければ、最初を
のように変更する。こうすると、左隅が0になり、下の隅が-3になる。幅は変わらないので、x座標は0から20に、y座標は-3から11ということになる。
基本的な図形を描くためには
ps.writePoint(x,y,r,col) | 点(●)を(x,y)に半径r、色colで描く |
ps.writeLine(x1,y1,x2,y2,col,w) | (x1,y1)から(x2,y2)に色colで幅wの線を引く。 |
ps.writeCircle(x,y,r,col,w) | (x,y)を中心として半径rの円を色colで描く。 |
ps.fillCircle(x,y,r,col,w) | 上に同じだが中身を塗りつぶす |
ps.writeRect(x1,y1,x2,y2,col,w) | (x1,y1)(x2,y2)を角にする長方形を描く。 |
ps.fillRect(x1,y1,x2,y2,col,w) | 上に同じだが中身を塗りつぶす |
ps.writeCoordinateXY() | x,y座標系を描く |
ps.writeCoordinateX() | x座標系のみ描く |
ps.writeCoordinateY() | y座標系のみ描く |
などの命令があるので利用してもよい。
javascriptでcanvasに図を描く一般的な方法については、Mozilla Developer Network(MDN)の「canvasに図形を描く」のページなどを参照のこと。ただし、このページで説明されている方法で使われている「グラフィックコンテキスト」は、ps.ctx(Phystemの名前がpsでないのなら、ps.の部分を書き換えること)で定義されている。
Phystemは「系」である。物理における「系」とは互いに相互作用するいくつかの物体(このPhystemで使えるのは質点のみ)の集まりであるが、最初にnew Phystemした段階では、物体は何もない。これから物体を系に登録していかなくてはいけない。
Phystemに「登録」することができる物体はDynamicalObjectとNonDynamicalObjectの2種類(と、これらを「継承」したオブジェクト)がある。DynamicalObjectは名前の通り「動力学的」な「もの」で、運動方程式に従って運動している物体である。NonDynamicalObjectの方は動いていないか、動いているとしても運動方程式によらない運動をしている物体である。Phystemによるシミュレーションの準備は、
と行われる。その後は物理法則にしたがって物体が運動する。
物体の登録は、
のように行う。最初のpsが登録される「系」である。x,yが物体の初期位置、vx,vyが物体の初速度であり、mは質量である。最後の"rgba(255,255,0,0.5)"は色を指定している。
系にDynamicalObjectが登録されたことにより、以後このDynamicalObjectは一定時間ごとにそれぞれの物理法則にしたがって『次の時刻での位置』を計算され、そうやって動いた結果がcanvasに表示される。
var m=new DynamicalObject(ps,x,y,vx,vy,m,"rgba(255,255,0,0.5)");を行うと、Phystemが持っているdyObjsという配列にこのmが登録される。
dyObjsに登録されていると、
という処理が行われる(1秒に何回か、決まった回数)。
.draw()や.drawV()などは物体のプロパティとして定義された関数で、変更することもできる(しなければm.rで設定された半径の円がm.colで設定された色で描かれる)。
Phystemがどのように力学的計算を行っているか、その部分を抜粋しよう。
for (i = 0; i < N; i++) { now = this.dyObjs[i]; F = new Vector(0, 0); for (j = 0 ; j < N; j++) { // j番目の動的オブジェクトからi番目の動的オブジェクトへの力 if (i !== j) { f = this.dyObjs[j].interactionForce(now); F.add(f); } } for (j = 0 ; j < this.ndObjs.length ; j++) { // j番目の【非】動的オブジェクトからi番目の動的オブジェクトへの力 f = this.ndObjs[j].interactionForce(now); F.add(f); } now.v.add(F.prod(dt / now.mass)); // v = v+(F(x+vΔt)/m)Δt now.pos.copyFrom(now.npos); }
この中のnow = this.dyObjs[i];のnowが今着目している物体で、次のループではまずDynamicalObjectの一つ一つからこのnowに働く力をFという変数に足している。
その計算をしている部分が f = this.dyObjs[j].interactionForce(now);の部分で、interactionForce()という、this.dyObjs[j](←j番目に登録されたDynamicalObject)のプロパティになっている関数を(nowを引数に)呼び出して力を計算させている。つまり、「this.dyObjs[j]がnowにどんな力を及ぼすか」はthis.dyObjs[j]というオブジェクトが知っている、というわけである。
この後次に、登録されているNonDynamicalObject(this.ndObjs[j])に対しても同様に、nowに対してどのような力を及ぼすかを計算する。これで系にある全ての物体からの力を計算したことになる。
ところで、運動方程式は
という式であるが、これを
のように連立微分方程式として解いている。具体的には
という手順を行っている。
先に位置を更新して、更新した位置を下に速度を更新するのは、「シンプレクティック数値積分法」という数値積分の方法で、オイラーの方法やルンゲ・クッタの方法に比べエネルギー保存則の成立において有利だと言われている。
この二つはphystem.jsの中で定義されている。phystem.jsの最初の方を見るとまず、
var PhysObject = function (ps, x, y, c, rr) { this.pos = new Vector( (x === undefined) ? 0 : x, (y === undefined) ? 0 : y ); // 色をセット this.col = (c === undefined) ? "rgba(0,0,0,0.4)" : c; this.npos = this.pos.makeCopy(); // npは、「次の位置」(シンプレクティック法で微分方程式を解いているのでこれが必要) this.r = (rr === undefined) ? 0.5 :rr; this.ps = ps; // 自分の属している系をセットしておく。 // この「系」は描画用のキャンバスも含んでいる。 }; PhysObject.prototype = { // 描画:●を描く。 draw: function () { this.ps.fillCircle(this.pos.x, this.pos.y, this.r, this.col); }, // 相互作用力。「to」へと及ぼす力を返す。 // 現在はスタブなので0ベクトルを返す。 interactionForce: function (to) { return new Vector(0, 0); }, // 変位べクトル。現在位置ではなく、△t後の位置で計算する。 displacementFrom: function (from) { return this.npos.diff(from.npos); } };
という部分がある。PhysObjectはDynamicalObjectとNonDynamicalObjectに共通する属性を持ったクラスで、このあとPhysObjectはDynamicalObjectとNonDynamicalObjectに「枝分かれ」する。
PhysObject自体は、this.pos(現在位置のベクトル)、this.npos(しばらく後での位置のベクトル)なぜこれがいるのかは後で説明する。、this.col(色)、this.ps(自分がどの系に属しているか)、this.r(自分の半径:円だと考えている)だけを属性(プロパティ)として持つ。
NonDynamicalObjectのコンストラクタは、
// 非力学オブジェクト // 運動しないもしくは運動方程式によらない運動をするオブジェクト var NonDynamicalObject = function (ps, x, y, c, r) { PhysObject.call(this, ps, x, y, c, r); if (ps !== undefined) { ps.ndObjs.push(this); // 「系」に、自分自身を登録する。 } }; NonDynamicalObject.prototype = Object.create(PhysObject.prototype); NonDynamicalObject.prototype.constructor = NonDynamicalObject;
のようになっている(最後の2行はjavascriptでクラスの継承を行うときのおまじないで、これでPhysObjectの持っているプロトタイプを受け継ぐことができる)。
ps.ndOjbs.push(this);は系(ps)の持っているndOjbsという名前のリストに自分自身(this)を登録している。これで、上で説明したループの中でこの物体が他の物体に及ぼす力が計算されるようになる。
一方、DynamicalObjectの方は、
// DynamicalObject:「動く物体」のクラス // 親クラスはPhysObject(物理的物体?) // mass(質量)と、fList(働いている力のリスト)という属性が追加される。 // yavというオブジェクトを持っているが、外部から使う必要は多分ない。 // 初期位置、初速度、質量、色、と「どの系に属すか」が呼び出しパラメータ(初速度以降は省略可) var DynamicalObject = function (ps, x, y, vx, vy, m, c, r) { PhysObject.call(this, ps, x, y, c, r); if (ps !== undefined) { // 自分自身を系の「動くものリスト」に登録する。 ps.dyObjs.push(this); } this.v = new Vector( (vx === undefined) ? 0 : vx, (vy === undefined) ? 0 : vy ); this.mass = (m === undefined) ? 1 : m; this.makeV(); // 動くものは速度や力を表示するので、その為の矢印を先に作っておく。 this.fList = []; // 力のリスト。 };
のようにコンストラクタが定義されている(おまじない部分もあるが、省略した)。
ps.dyOjbs.push(this);で、今度は系の持っているdyObjsに自分自身を登録する。Dynamicalなので初速度や、働く力など、他の量も最初に用意している。
上で述べたように、.interactionForce()という関数が「この物体が他の物体にどんな力を及ぼすか」を設定しているので、ここを書き換えることで物理法則を変更できる。
最初の段階では、phystem.jsの中の
interactionForce: function (to) { return new Vector(0, 0); },
の部分でPhysObjectが出す力が「0ベクトル」に設定されている。
「書き換える」と言ってもここを変えてしまうと、全ての物体の力がいっせいに変わってしまう。たとえばある物体mの「力」を変更するには、
m.interactionForce=function(to) { mからtoに働く力を計算する部分 return 結果のVector; }
のようなコードを付け足せばよい。
しかし、これだとたくさんの物体に力を設定するのがたいへんである。そこでここでもオブジェクト指向プログラミングの恩恵を使うことにする。
実例で示そう。
表紙でも見せていた分子運動シミュレーションでは、
var Molecule=function (ps,x,y,vx,vy,m,c) { DynamicalObject.call(this,ps,x,y,vx,vy,m,c); }; Molecule.prototype=Object.create(DynamicalObject.prototype); Molecule.prototype.constructor=Molecule; Molecule.prototype.interactionForce=function (to) { var r=to.npos.diff(this.npos); var rlen=r.length(); return r.prod(10*Math.exp(-2*rlen*rlen)); };
というコードが書かれている。
1行目のvar Molecule=function (ps,x,y,vx,vy,m,c) {から始まる部分は、Moleculeという新しいオブジェクトが作られた時にどのような関数が呼ばれるかが書かれている。その中身はDynamicalObject.call(this,ps,x,y,vx,vy,m,c);で、これは「DynamicalObjectが作られた時に呼ばれる関数を実行してくれ」という、言わば「丸投げ」である。
その後に書かれている、
Molecule.prototype=Object.create(DynamicalObject.prototype); Molecule.prototype.constructor=Molecule;
の部分は、MoleculeのプロトタイプをDynamicalObjectのプロトタイプと一致させるための命令である(詳細は省くので呪文と思って実行すること)。
その後のMolecule.prototype.interactionForce=function (m2) {で始まる部分がいよいよ力を定義する部分で、Molecule.protypeとあることで「Moleculeのプロトタイプ」にinteractionForceを登録する(プロトタイプに登録することで、以後new Moleculeで作られたオブジェクトにはこのinteractionForceが登録された状態になる)。
力の中身についても説明しておこう。
var r=to.npos.diff(this.npos);では、rという変数にtoのnpos(次の段階での位置)とthis(自分自身)のnposとの差のベクトルを設定している。
次のvar rlen=r.length();はそのベクトルの長さをrlenという新しい変数に代入した。
r.prod(10*Math.exp(-2*rlen*rlen))は、差のベクトルrを(10*Math.exp(-2*rlen*rlen)倍して、これを働く力として返している。こうすることで自分から相手に向かう方向(図の緑矢印)を向いた力(大きさは最後の掛算で調整している)をtoに対して与えることになる。
新しいタイプの「物体」を作りたかったら、以下の手順を実行する。
var 新クラス=function (ps,x,y,vx,vy,m,c) { DynamicalObject.call(this,ps,x,y,vx,vy,m,c); }; 新クラス.prototype=Object.create(DynamicalObject.prototype); 新クラス.prototype.constructor=新クラス;を書く。新クラスの部分は自分で決めた新しいクラスの名前を入れること。
以下は、万有引力を及ぼし合う惑星を定義する場合のinteractionForceの例である。
Planet.prototype.interactionForce=function (m2) { var r=m2.npos.diff(this.npos); var rlen=r.length(); var rdis=m2.r + this.r; if( rlen > rdis ) { return r.prod(-G*this.mass*m2.mass/ (rlen*rlen*rlen) ); } else { return r.prod(-G*this.mass*m2.mass/ (rdis*rdis*rdis) ); } };
ここではDynamicalObjectの場合を考えたが、NonDynamicalObjectも(他の物体に力を及ぼすことはできるので)同様に力を定義できる。例として、ps.makeGravity();で重力を作っている方法を説明しよう。
makeGravityは以下のように定義されている(phystem.jsの中にある)。
// 重力を追加する Phystem.prototype.makeGravity = function (g) { this.gravity = (g === undefined) ? 1 : g; var earth = new NonDynamicalObject(this); earth.interactionForce = function (to) { return new Vector(0, -to.mass * to.ps.gravity); }; earth.draw = function () {}; };
ここでやっていることは「地球(earth)」というNonDynamicalObjectを作って(一個しか作ってないのでクラスを定義していない)、そのearthのinteractionForceを、Vector(0, -to.mass * to.ps.gravity)にしている(これは下向きmgを表している)。これでPhystemに「earth」という非力学的オブジェクトが登録され、earthからの力として重力が他の全ての物体にかかることになる。
なお最後のearth.draw = function () {};では、地球を描画する必要はないので、描画する関数を「何もしない{}」にセットしている。
makeAir();(←空気抵抗を追加する)も同様なので、phystem.jsを見てみよう。
前回で見せた2個の人工衛星の追っかけっ子の中では、
M=new NonDynamicalObject(gs,0,0,"rgba(0,200,0,0.8)"); M.interactionForce=function (to) { var r=to.npos; var rlen=r.length(); if (rlen > 0.5) { return r.quot(-(rlen*rlen*rlen)); } else { return r.prod(-8); } }
のようにして、地球を表すオブジェクトMが人工衛星に万有引力を及ぼすようにした。
他に、NonDynamicalObjectの例としてはBaneもある(バネは動くのだが、その動きは運動方程式に応じたものではないので、NonDynamicalである)。bane.jsを見てみるとよい。
// m1とm2をつなぐバネ var Bane = function (ps, m1, m2, k, l, c) { // バネはNonDynamicalObjectである。 NonDynamicalObject.call(this, ps, 0, 0, c); this.m1 = m1; this.m2 = m2; // バネ定数(デフォルトは1) this.k = (k === undefined) ? 1 : k; // 自然長(デフォルトは0) this.l = (l === undefined) ? 0 : l; };
まずBaneはNonDynamicalObjectの一種なので、NonDynamicalObjectのコンストラクタを呼び出す(この後に、プロトタイプを一致させるためのいつもの呪文もある)。
バネは二つのDynamicalObjectをつないでいるものなので、その二つをthis.m1とthis.m2に格納する。その後はバネ定数と自然長をセットしている。
Baneのプロトタイプがその後につづいているが、
Bane.prototype.interactionForce = function (to) { if (to === this.m1) { return this.force1(); } else if (to === this.m2) { return this.force2(); } else { return new Vector(0, 0); } };
の部分で.interactionForce()を定義しているが、これはto(力を呼ぼす相手)がm1かm2だった時にだけ0でない力を返すようになっている(つながっていない物体にまでバネが力を出したらたいへんだ)。力の計算自体は.force1()と.force2()という別の(やはりプロトタイプで定義された)関数に丸投げされている(どうなっているか知りたい人はソースを確認すること)。