用Processing练习写的一个椭圆运动模拟,代码为java格式,完整代码附在下文。
在setup函数中配置行星位置(init_pos)和初速度(init_vel)。设初速度为(x, y),当y = 0时可以模拟圆周运动,其他情况可以模拟椭圆运动。
当椭圆的离心率很大时,会出现下图的轨道进动的现象,这实际上是由误差累积导致的:模拟程序初始化时是设置的偏远日点的速度,椭圆轨道近日点速度和远日点速度约成反比,因此离心率越高,近日点速度越大。而软件计算和更新位置是恒定速度的,速度过快导致微小计算误差通过乘以大的速度放大,从而导致了轨道偏移。这也给了一条经验:轨道模拟时需要防止近日点速度过大,或者提高计算精度。
int MASS_SIZE_RATE = 2;
class Mover {
PVector pos;
PVector last_pos;
PVector vel;
PVector acc = new PVector(0, 0);
float mass;
Mover(PVector init_pos, PVector init_vel, float mass) {
pos = init_pos;
last_pos = init_pos.copy();
vel = init_vel;
this.mass = mass;
}
void applyForce(PVector force) {
acc.add(force.div(mass));
}
void update() {
vel.add(acc);
pos.add(vel);
acc.x = 0;
acc.y = 0;
}
void display() {
stroke(0);
// fill(175);
// println(last_pos.x, last_pos.y, pos.x, pos.y);
// ellipse(pos.x, pos.y, mass * MASS_SIZE_RATE, mass * MASS_SIZE_RATE);
//float line_vel = vel.mag();
//println(line_vel);
line(last_pos.x, last_pos.y, pos.x, pos.y);
last_pos.x = pos.x;
last_pos.y = pos.y;
}
}
class Attractor {
PVector pos;
float mass;
float G;
Attractor(PVector init_pos, float mass, float G) {
pos = init_pos;
this.mass = mass;
this.G = G;
}
void display() {
stroke(0);
fill(255);
ellipse(pos.x, pos.y, mass * MASS_SIZE_RATE, mass * MASS_SIZE_RATE);
}
PVector attract(Mover m) {
PVector r = PVector.sub(pos, m.pos);
float distance = r.mag();
PVector direction = r.normalize();
float f = (G * mass * m.mass) / (distance * distance);
return direction.mult(f);
}
}
Mover mover;
Attractor attr;
void setup() {
fullScreen();
background(255);
PVector init_pos = new PVector(width/2 - 400, height/2);
PVector init_vel = new PVector(1.5, 1);
init_vel.mult(2);
mover = new Mover(init_pos, init_vel, 1);
PVector attr_pos = new PVector(width / 2, height / 2);
attr = new Attractor(attr_pos, 10, 1000);
}
void draw() {
PVector f = attr.attract(mover);
mover.applyForce(f);
mover.update();
mover.display();
attr.display();
}