Processing 椭圆运动模拟

用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();
}