惑星シミュレーション

[いちごライス]

下記のコードは現在作成中の惑星シミュレーションのコードです。このコードに粒子間の分子間力を組み込みたいのですが、何をどうすれば良いか分からず、手詰まっております。物理に詳しいかたがもしいらっしゃったら、手を貸して頂きたいです。

なお、作ろうとしているのは、真ん中に重い星を置き、その周りを多数の粒子からなる惑星が公転運動している。というものです。




#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

//定義
#define PI 3.1415926535
#define GM 4 * PI *PI
#define M 1.0 // 惑星の質量
#define dT 1.0 / 256.0 // 時間ステップ
#define X0 1.0 // x 方向の位置の初期値
#define Y0 0.0 // y 方向の位置の初期値
#define VX0 0.0 // x 方向の速度の初期値
#define VY0 sqrt(GM) * M // y 方向の速度の初期値
#define RADIUS 0.05 // 惑星の半径
#define INIROLL 0.0
#define ENDROLL 1.0
#define NA 3 //粒子の数

//クラスの作成
class Planet {
public:
double x, y;
double vx, vy;

void M_IC(); //Main & Initial Condition
void RK(); //Runge-Kutta
void print();

};

// x,yの初期値
void Planet::M_IC() {

double xa, ya;
double sq;


srand((unsigned)time(NULL)); // 乱数の仕組みを初期化

for (int i = 0; i < NA;) {
printf("------ R[%d] ------\n", i);
xa = (double)rand() / RAND_MAX * 2 * RADIUS - RADIUS;
ya = (double)rand() / RAND_MAX * 2 * RADIUS - RADIUS;
sq = xa * xa + ya * ya; //半径の２乗

if (sq > RADIUS) {
continue;
} // sq が半径を超えている点は不採用
i++;

x = xa + X0;
y = ya + Y0;

vx = VX0;
vy = VY0;


RK();
print();
}

}


//ルンゲクッタ＆csvファイルへの書き込み
void Planet::RK() {

FILE *fp;
double xk1, yk1, vxk1, vyk1;
double xk2, yk2, vxk2, vyk2;
double xk3, yk3, vxk3, vyk3;
double xk4, yk4, vxk4, vyk4;
double tx, ty, tvx, tvy;
double q, qi3;

fopen_s(&fp, "TD.csv", "a");// csvファイルへの書き込み

for (double t = INIROLL; t < ENDROLL; t += dT) {

tx = x;
ty = y;
tvx = vx;
tvy = vy;
q = hypot(tx, ty);
qi3 = 1.0 / (q * q * q);
xk1 = tvx / M * dT;
yk1 = tvy / M * dT;
vxk1 = -GM * M * tx * qi3 * dT;
vyk1 = -GM * M * ty * qi3 * dT;

tx = x + 0.5 * xk1;
ty = y + 0.5 * yk1;
tvx = vx + 0.5 * vxk1;
tvy = vy + 0.5 * vyk1;
q = hypot(tx, ty);
qi3 = 1.0 / (q * q * q);
xk2 = tvx / M * dT;
yk2 = tvy / M * dT;
vxk2 = -GM * M * tx * qi3 * dT;
vyk2 = -GM * M * ty * qi3 * dT;

tx = x + 0.5 * xk2;
ty = y + 0.5 * yk2;
tvx = vx + 0.5 * vxk2;
tvy = vy + 0.5 * vyk2;
q = hypot(tx, ty);
qi3 = 1.0 / (q * q * q);
xk3 = tvx / M * dT;
yk3 = tvy / M * dT;
vxk3 = -GM * M * tx * qi3 * dT;
vyk3 = -GM * M * ty * qi3 * dT;

tx = x + xk3;
ty = y + yk3;
tvx = vx + vxk3;
tvy = vy + vyk3;
q = hypot(tx, ty);
qi3 = 1.0 / (q * q * q);
xk4 = tvx / M * dT;
yk4 = tvy / M * dT;
vxk4 = -GM * M * tx * qi3 * dT;
vyk4 = -GM * M * ty * qi3 * dT;

x += (xk1 + 2 * xk2 + 2 * xk3 + xk4) * (1.0 / 6);
y += (yk1 + 2 * yk2 + 2 * yk3 + yk4) * (1.0 / 6);
vx += (vxk1 + 2 * vxk2 + 2 * vxk3 + vxk4) * (1.0 / 6);
vy += (vyk1 + 2 * vyk2 + 2 * vyk3 + vyk4) * (1.0 / 6);

print();
fprintf(fp, "%f, %f\n", x, y);

}
fclose(fp);
}

//出力
void Planet::print() {
printf("x=%f y=%f vx=%f vy=%f \n", x, y, vx, vy);
}

//メインプログラム
int main() {
Planet R[NA];

R[NA].M_IC();


return 0;
}

[幸尚]

こんばんは、物理のことは全然分かりませんがソースコード上に「ルンゲクッタ」というキーワードがあったので
ググってみたところたくさんの計算方法などがヒットしました。
まず、ルンゲクッタというのを調べてみてはどうでしょうか？
なにか手がかりが見つかるかもしれませんよ。
また、「ルンゲクッタ法 C言語」でググってみたところ同様にたくさんのサンプルなどがありました。

[usao]

示されたコードに関する詳細な状況は不明ですが，
ぱっと見の雰囲気的には（ここまで書いておいていまさら）「ルンゲクッタ法とは何か？」とかいうことではないであろう，と想像します．

何が不明点なのかも不明ですが，

> 多数の粒子からなる惑星

っていうのを「少数の粒子で」まずはやってみればよいのでは．
いきなり多数だと，挙動の確認も無理（：「こうなるはず」っていう相互作用が人間に把握困難になる）でしょうし．

何を持って「分子間力」とするのか謎ですが，引力と斥力みたいなのを相応にモデル化できれば，あとはそれによる運動を差分法で扱ってやるだけなのでは．

[あたっしゅ]

東上☆海美☆「
#1＞ 下記のコードは現在作成中の惑星シミュレーションのコードです。
#1＞ このコードに粒子間の分子間力を組み込みたいのですが、何をどうすれば良いか分からず、手詰まっております。

https://ja.wikipedia.org/wiki/%E5%88%86 ... 3%E5%8A%9B
分子間力 - ジャパリ図書館(ja)

https://ja.wikipedia.org/wiki/%E5%BC%95 ... 5%E5%8A%9B
引力と斥力 - ジャパリ図書館(ja)

『粒子の分子間力』と『惑星を含む天体間の引力』とは、全く異なるものみみ。


#1＞ なお、作ろうとしているのは、真ん中に重い星を置き、その周りを多数の粒子からなる惑星が公転運動している。というものです。

貴方のプログラムは、3 つの惑星間の関係について計算しようとしているようです。
『真ん中に重い星(太陽)』は、無いし、『多数の粒子からなる惑星』は、実装されていません。

もしかして、多数の粒子が集まって、一つの惑星になるシミュレーションをしたいのですか ?

普通、『惑星の公転シミュレーション』と言うと、『真ん中に重い星(太陽)』と惑星ひとつひとつについて
注目して、惑星間の引力は無視しますが。
」