magnus.cpp

custom_math::vector_3 acceleration(custom_math::vector_3 pos, custom_math::vector_3 vel, custom_math::vector_3 ang_vel)
{
	// http://twu.tennis-warehouse.com/learning_center/aerodynamics2.php
	const double air_density = 1.225;
	const double lift_coeff = 0.05;
	const double drag_coeff = 0.55;
	const double ball_cross_section = 0.0034;
	const double ball_mass = 0.0585;

	// Gravitation, in metres per second, per second
	custom_math::vector_3 grav_accel(0, -9.81, 0);

	// Magnus effect, in metres per second, per second
	// http://farside.ph.utexas.edu/teaching/329/lectures/node43.html
	// http://spiff.rit.edu/richmond/baseball/traj/traj.html
	custom_math::vector_3 magnus_accel = vel.cross(ang_vel)*0.5*air_density*lift_coeff*ball_cross_section / ball_mass;

	// Wind and drag, in metres per second, per second
	custom_math::vector_3 wind_vel(2.77778, 0, 0); // in metres per second
	custom_math::vector_3 drag_vel = wind_vel - vel;
	custom_math::vector_3 drag_accel = drag_vel*drag_vel.length()*0.5*air_density*drag_coeff*ball_cross_section / ball_mass;
	
	return grav_accel + magnus_accel + drag_accel;
}

// https://en.wikipedia.org/wiki/Symplectic_integrator
// Also known as Verlet integration
void proceed_symplectic2(custom_math::vector_3 &pos, custom_math::vector_3 &vel, const custom_math::vector_3 &ang_vel)
{
	static const double c[2] = { 0, 1 };
	static const double d[2] = { 0.5, 0.5 };

// 	pos += vel * c[0] * dt; // first element c[0] is always 0
	vel += acceleration(pos, vel, ang_vel) * d[0] * dt;

	pos += vel * c[1] * dt;
	vel += acceleration(pos, vel, ang_vel) * d[1] * dt;
}

// https://www.gamedev.net/forums/topic/701376-weird-circular-orbit-problem/?do=findComment&comment=5402054
// https://en.wikipedia.org/wiki/Symplectic_integrator
void proceed_symplectic4(custom_math::vector_3 &pos, custom_math::vector_3 &vel, const custom_math::vector_3 &ang_vel)
{
	static double const cr2 = pow(2.0, 1.0 / 3.0);

	static const double c[4] =
	{
		1.0 / (2.0*(2.0 - cr2)),
		(1.0 - cr2) / (2.0*(2.0 - cr2)),
		(1.0 - cr2) / (2.0*(2.0 - cr2)),
		1.0 / (2.0*(2.0 - cr2))
	};

	static const double d[4] =
	{
		1.0 / (2.0 - cr2),
		-cr2 / (2.0 - cr2),
		1.0 / (2.0 - cr2),
		0.0
	};

	pos += vel * c[0] * dt;
	vel += acceleration(pos, vel, ang_vel) * d[0] * dt;

	pos += vel * c[1] * dt;
	vel += acceleration(pos, vel, ang_vel) * d[1] * dt;

	pos += vel * c[2] * dt;
	vel += acceleration(pos, vel, ang_vel) * d[2] * dt;

	pos += vel * c[3] * dt;
//	vel += acceleration(pos, vel, ang_vel) * d[3] * dt; // last element d[3] is always 0
}

void proceed_Euler(custom_math::vector_3 &pos, custom_math::vector_3 &vel, const custom_math::vector_3 &ang_vel)
{
	vel += acceleration(pos, vel, ang_vel) * dt;
	pos += vel * dt;
}

inline void proceed_RK2(custom_math::vector_3 &pos, custom_math::vector_3 &vel, const custom_math::vector_3 &ang_vel)
{
	custom_math::vector_3 k1_velocity = vel;
	custom_math::vector_3 k1_acceleration = acceleration(pos, k1_velocity, ang_vel);
	custom_math::vector_3 k2_velocity = vel + k1_acceleration * dt*0.5;
	custom_math::vector_3 k2_acceleration = acceleration(pos + k1_velocity * dt*0.5, k2_velocity, ang_vel);

	vel += k2_acceleration * dt;
	pos += k2_velocity * dt;
}

void proceed_RK4(custom_math::vector_3 &pos, custom_math::vector_3 &vel, const custom_math::vector_3 &ang_vel)
{
	static const double one_sixth = 1.0 / 6.0;

	custom_math::vector_3 k1_velocity = vel;
	custom_math::vector_3 k1_acceleration = acceleration(pos, k1_velocity, ang_vel);
	custom_math::vector_3 k2_velocity = vel + k1_acceleration * dt*0.5;
	custom_math::vector_3 k2_acceleration = acceleration(pos + k1_velocity * dt*0.5, k2_velocity, ang_vel);
	custom_math::vector_3 k3_velocity = vel + k2_acceleration * dt*0.5;
	custom_math::vector_3 k3_acceleration = acceleration(pos + k2_velocity * dt*0.5, k3_velocity, ang_vel);
	custom_math::vector_3 k4_velocity = vel + k3_acceleration * dt;
	custom_math::vector_3 k4_acceleration = acceleration(pos + k3_velocity * dt, k4_velocity, ang_vel);

	vel += (k1_acceleration + (k2_acceleration + k3_acceleration)*2.0 + k4_acceleration)*one_sixth*dt;
	pos += (k1_velocity + (k2_velocity + k3_velocity)*2.0 + k4_velocity)*one_sixth*dt;
}