% -------------------------------------------------------------------------
% plot the path of a ball bouncing in a box
%
% physics in MKS units
% -------------------------------------------------------------------------

% video parameters
FileName = 'bounce.avi';
FrameRate = 10; % frames per second
KeyRate = 10; % key frames per second (=FrameRate => no interpolation)
Duration = 10; % total duration of animation in seconds

% simulation parameters
dt = .01; % time step per interval
ni = Duration / dt; % total number of intervals to calculate
np = 1 / FrameRate / dt; % number of intervals between frames


% physical parameters
rmin = [0 0 0]; rmax = [100 100 100]; % size of the box
rball = 10; % radius of the ball
m = 1; % mass of the ball
ke = 50; % constant of elasticity
ks = 2*ke/m*rball*dt; % energy stored in deformation
ks = repmat(ks,1,3); % make a vector out of it for calculation
kb = .015; % loss coefficient: fraction of energy lost per bounce
kd = sqrt(1 - kb) - 1; % velocity change at each bounce

vo = 10; alpha = pi/3; % initial horizontal velocity

r = [rball           rball           rmax(3)-rball]; % position
v = [vo*cos(alpha)   vo*sin(alpha)   30            ]; % velocity
a = [0               0              -9.81         ]; % acceleration

% plotting parameters
whitebg([0.97 0.99 1.00]); grid on;
ViewLimits = reshape([rmin; rmax], 1, 6);

% create the avi file and set its parameters
video = avifile( ...
  FileName,'fps', FrameRate, 'keyframe', KeyRate, ...
  'quality', 100, 'compression', 'None' ...
);

for i = 1:ni
  % detect collisions with the walls of the box
  r1 = r - rmin; c1 = (r1 < rball);
  r2 = rmax - r; c2 = (r2 < rball);
  % calculate radii of the possibly deformed shape
  rd = rball*(1-c1).*(1-c2) + r1.*c1 + r2.*c2;
  % update position and velocity
  r = r + v*dt;
  v = v + a*dt + c1.*(kd*v + ks./r1) + c2.*(kd*v - ks./r2);
  % plot position once every np intervals
  if rem(i, np) == 0
    % create the possibly deformed ball
    [Ex Ey Ez] = ellipsoid( ...
      r(1), r(2), r(3), rd(1), rd(2), rd(3) ...
    );
    % draw it and use the returned handle to set material properties
    h = surf(Ex, Ey, Ez);
    set(h, ...
      'FaceLighting','phong', 'FaceColor', [1,.9,.9], ...
      'EdgeLighting','phong', 'EdgeColor', 'none', ...
      'AmbientStrength',0.6, 'DiffuseStrength', 0.4, ...
      'SpecularStrength', 0.2, 'SpecularExponent', 8, ...
      'SpecularColorReflectance', 1.0 ...
    );
    % set lighting and view parameters
    light('Position', [1 0 1], 'Style', 'infinite');
    view(50,33); axis(ViewLimits,'square');
    % capture the current frame
    video = addframe(video, getframe);
  end
end
video = close(video);
