#1
by いちごライス » 1年前
下記のコードは現在作成中の惑星シミュレーションのコードです。このコードに粒子間の分子間力を組み込みたいのですが、何をどうすれば良いか分からず、手詰まっております。物理に詳しいかたがもしいらっしゃったら、手を貸して頂きたいです。
なお、作ろうとしているのは、真ん中に重い星を置き、その周りを多数の粒子からなる惑星が公転運動している。というものです。
#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; //半径の2乗
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;
}
下記のコードは現在作成中の惑星シミュレーションのコードです。このコードに粒子間の分子間力を組み込みたいのですが、何をどうすれば良いか分からず、手詰まっております。物理に詳しいかたがもしいらっしゃったら、手を貸して頂きたいです。
なお、作ろうとしているのは、真ん中に重い星を置き、その周りを多数の粒子からなる惑星が公転運動している。というものです。
#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; //半径の2乗
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;
}