anchor.diff:  % -------------------------------------------------------------------------
   
   -include("log.hrl").
  -%-include("pre_physics.hrl").
  -
  --record(physical, {
  -	% Updated values (assume these change every frame)
  -	position = {0, 0, 0} :: vector:vec(),
  -	linear_momentum = {0, 0, 0} :: vector:vec(),
  -	orientation = {1, 0, 0, 0} :: quaternion:quat(),
  -	angular_momentum = {0, 0, 0} :: vector:vec(),
  -
  -	% Input-only values
  -	force_abs = {0, 0, 0} :: vector:vec(),
  -	force_rel = {0, 0, 0} :: vector:vec(),
  -	torque_absolute = {0, 0, 0} :: vector:vec(),
  -	torque_relative = {0, 0, 0} :: vector:vec(),
  -
  -	% Purely calculated values (DON'T try to change these externally)
  -	last_update :: erlang:timestamp(),
  -	linear_velocity = {0, 0, 0} :: vector:vec(),
  -	angular_velocity = {0, 0, 0} :: vector:vec(),
  -	spin = {1, 0, 0, 0} :: quaternion:quat(),
  -
  -	% Constant values
  -	mass = 1 :: float(),
  -	inverse_mass = 1 :: float(),
  -	inertia_tensor = 1 :: float(),
  -	inverse_inertia_tensor = 1 :: float()
  -}).
  +-include("pre_physics.hrl").
   
   %% ------------------------------------------------------------------------
   %% External API
  
positive.code: 	#physical{
  		position = Position,
  		linear_velocity = Velocity,
  		orientation = Orientation,
  		angular_momentum = AngularMomentum
  	} = State,
  	
  	{Velocity1, Force1, Spin1, Torque1, State1} = evaluate(0, 0, 0, 0, 0, State),
  	{Velocity2, Force2, Spin2, Torque2, State2} = evaluate(TimeDelta * 0.5, Velocity1, Force1, Spin1, Torque1, State1),
  	{Velocity3, Force3, Spin3, Torque3, State3} = evaluate(TimeDelta * 0.5, Velocity2, Force2, Spin2, Torque2, State2),
  	{Velocity4, Force4, Spin4, Torque4, State4} = evaluate(TimeDelta, Velocity3, Force3, Spin3, Torque3, State3),
  
  	NewVelocity = 1.0 / 6.0 * (Velocity1 + 2.0 * (Velocity2 + Velocity3) + Velocity4),
  	NewAcceleration = 1.0 / 6.0 * (Force1 + 2.0 * (Force2 + Force3) + Force4),
  	NewSpin = 1.0 / 6.0 * (Spin1 + 2.0 * (Spin2 + Spin3) + Spin4),
  	NewTorque = 1.0 / 6.0 * (Torque1 + 2.0 * (Torque2 + Torque3) + Torque4),
  
  	State4#physical{
  		position = Position + NewVelocity * TimeDelta,
  		linear_velocity = Velocity + NewAcceleration * TimeDelta,
  		orientation = quaternion:unit(Orientation + NewSpin * TimeDelta),
  		angular_momentum = AngularMomentum + NewTorque * TimeDelta
  	}.
  
negative.code: 		angular_momentum = AngularMomentum
  	} = State,
  
  	{SpinW, SpinX, SpinY, SpinZ} = Spin,
  	{OrientW, OrientX, OrientY, OrientZ} = Orientation,
  	NextPosition = vector:add(Position, vector:multiply(TimeDelta, Velocity)),
  	NextVelocity = vector:add(InitialVelocity, vector:multiply(TimeDelta, Force)),
  	%XXX: I have no idea if/how this works, but it's what the example code did...
  	NextOrientation = {
  		OrientW + SpinW * TimeDelta,
  		OrientX + SpinX * TimeDelta,
  		OrientY + SpinY * TimeDelta,
  		OrientZ + SpinZ * TimeDelta
  	},
  	NextAngularMomentum = vector:add(AngularMomentum, vector:multiply(TimeDelta, Torque)),
  
  	State1 = State#physical{
  		position = NextPosition,
  		linear_velocity = NextVelocity,
  		orientation = quaternion:unit(NextOrientation),
  		angular_momentum = NextAngularMomentum
  	},
  	%TODO: Should _NextAngularVelocity even be returned?
  	{NextSpin, State2} = update_state(State1),
  
  	{NextForce, NextTorque} = forces(TimeDelta, State2),
  
  	{NextVelocity, NextForce, NextSpin, NextTorque, State2}.
  
  %% ------------------------------------------------------------------------
  
  