cjameshuff/gist:3909805

## gistfile1.cpp
struct Body {
    body_id_t id;

    double mass;

    vec2 pos, tmp_pos;
    vec2 vel, tmp_vel;
    vec2 accel, prev_accel;
    vec2 kv1, kv2, kv3, kv4;
    vec2 kp1, kp2, kp3, kp4;

    Body() {}

    Body(body_id_t bid, double m, const vec2 & p, const vec2 & v):
        id(bid),
        mass(m),
        pos(p),
        vel(v)
    {}

    // Copy constructor: only copies primary data, not intermediate values
    Body(const Body & rhs):
        id(rhs.id),
        mass(rhs.mass),
        pos(rhs.pos),
        vel(rhs.vel)
    {}
};

void System::StepRK4(double tstep)
{
    const double tstep_2 = tstep/2.0;
    const double tstep_6 = tstep/6.0;
    std::vector<Body *>::iterator b;
    // Save initial position and velocity
    for(b = bodies.begin(); b != bodies.end(); ++b) {
        (*b)->tmp_pos = (*b)->pos;
        (*b)->tmp_vel = (*b)->vel;
    }

    // k1 estimate
    DoForces();
    for(b = bodies.begin(); b != bodies.end(); ++b) {
        (*b)->kv1 = (*b)->accel;
        (*b)->kp1 = (*b)->vel;
    }

    // k2 estimate
    for(b = bodies.begin(); b != bodies.end(); ++b)
        (*b)->pos = (*b)->tmp_pos + (*b)->kp1*tstep_2;
    DoForces();
    for(b = bodies.begin(); b != bodies.end(); ++b) {
        (*b)->kv2 = (*b)->accel;
        (*b)->kp2 = (*b)->tmp_vel + (*b)->kv1*tstep_2;
    }

    // k3 estimate
    for(b = bodies.begin(); b != bodies.end(); ++b)
        (*b)->pos = (*b)->tmp_pos + (*b)->kp2*tstep_2;
    DoForces();
    for(b = bodies.begin(); b != bodies.end(); ++b) {
        (*b)->kv3 = (*b)->accel;
        (*b)->kp3 = (*b)->tmp_vel + (*b)->kv2*tstep_2;
    }

    // k4 estimate
    for(b = bodies.begin(); b != bodies.end(); ++b)
        (*b)->pos = (*b)->tmp_pos + (*b)->kp3*tstep;
    DoForces();
    for(b = bodies.begin(); b != bodies.end(); ++b) {
        (*b)->kv4 = (*b)->accel;
        (*b)->kp4 = (*b)->tmp_vel + (*b)->kv3*tstep;
    }

    // Update position and velocity
    for(b = bodies.begin(); b != bodies.end(); ++b) {
        vec2 ksum = (*b)->kv1 + 2.0*((*b)->kv2 + ((*b)->kv3)) + (*b)->kv4;
        (*b)->vel = (*b)->tmp_vel + ksum*tstep_6;

        ksum = (*b)->kp1 + 2.0*((*b)->kp2 + ((*b)->kp3)) + (*b)->kp4;
        (*b)->pos = (*b)->tmp_pos + ksum*tstep_6;
    }
}
	struct Body {
	body_id_t id;

	double mass;

	vec2 pos, tmp_pos;
	vec2 vel, tmp_vel;
	vec2 accel, prev_accel;
	vec2 kv1, kv2, kv3, kv4;
	vec2 kp1, kp2, kp3, kp4;

	Body() {}

	Body(body_id_t bid, double m, const vec2 & p, const vec2 & v):
	id(bid),
	mass(m),
	pos(p),
	vel(v)
	{}

	// Copy constructor: only copies primary data, not intermediate values
	Body(const Body & rhs):
	id(rhs.id),
	mass(rhs.mass),
	pos(rhs.pos),
	vel(rhs.vel)
	{}
	};

	void System::StepRK4(double tstep)
	{
	const double tstep_2 = tstep/2.0;
	const double tstep_6 = tstep/6.0;
	std::vector<Body *>::iterator b;
	// Save initial position and velocity
	for(b = bodies.begin(); b != bodies.end(); ++b) {
	(b)->tmp_pos = (b)->pos;
	(b)->tmp_vel = (b)->vel;
	}

	// k1 estimate
	DoForces();
	for(b = bodies.begin(); b != bodies.end(); ++b) {
	(b)->kv1 = (b)->accel;
	(b)->kp1 = (b)->vel;
	}

	// k2 estimate
	for(b = bodies.begin(); b != bodies.end(); ++b)
	(b)->pos = (b)->tmp_pos + (b)->kp1tstep_2;
	DoForces();
	for(b = bodies.begin(); b != bodies.end(); ++b) {
	(b)->kv2 = (b)->accel;
	(b)->kp2 = (b)->tmp_vel + (b)->kv1tstep_2;
	}

	// k3 estimate
	for(b = bodies.begin(); b != bodies.end(); ++b)
	(b)->pos = (b)->tmp_pos + (b)->kp2tstep_2;
	DoForces();
	for(b = bodies.begin(); b != bodies.end(); ++b) {
	(b)->kv3 = (b)->accel;
	(b)->kp3 = (b)->tmp_vel + (b)->kv2tstep_2;
	}

	// k4 estimate
	for(b = bodies.begin(); b != bodies.end(); ++b)
	(b)->pos = (b)->tmp_pos + (b)->kp3tstep;
	DoForces();
	for(b = bodies.begin(); b != bodies.end(); ++b) {
	(b)->kv4 = (b)->accel;
	(b)->kp4 = (b)->tmp_vel + (b)->kv3tstep;
	}

	// Update position and velocity
	for(b = bodies.begin(); b != bodies.end(); ++b) {
	vec2 ksum = (b)->kv1 + 2.0((b)->kv2 + ((b)->kv3)) + (*b)->kv4;
	(b)->vel = (b)->tmp_vel + ksum*tstep_6;

	ksum = (b)->kp1 + 2.0((b)->kp2 + ((b)->kp3)) + (*b)->kp4;
	(b)->pos = (b)->tmp_pos + ksum*tstep_6;
	}
	}