% bounce_credits Demonstration script for HW3.4 [PHYS210-HW3]
% 
% Full-blown version of bounce including recording of 
% cheesy credits.
% 
%***********************************************************
% IMPORTANT!!
%
%    DO NOT base mbounce.m on this file.  Start from bounce.m
%    instead.
% 
%    The TAs will NOT want deal with the volume of comments 
%    contained in this version.
%
%    Note that bounce.m doesn't add credits to the AVIs by 
%    so if you want to add those, refer to the last section 
%    of the code to see my approach, which you are free
%    to adopt.
%***********************************************************
%
% This script implements the "bouncing ball" simulation 
% described in the handout.
%  
% Its primary purpose is to serve as a template for your 
% solution of Problem 4.
% 
% Its secondary purpose is to illustrate the following:
%  
% (1) A simple Matlab animation technique wherein at each
%     time step of the simulation the plot window is cleared
%     and all elements of the plot are redrawn.
%  
% (2) The creation of an AVI format video file that records
%     the animation.
%  
% (3) Some additional features of the plot command that aren't
%     discussed in the lab notes or course script/functions 
%     files.
%  
%-----------------------------------------------------------
% DESCRIPTION OF SIMULATION
%-----------------------------------------------------------
% 
% Refer to the handout for the full description.
%  
% The ball (particle) is launched with some initial 
% velocity and is accelerated downwards by gravity above a 
% floor with infinitesimal thickness in the y-direction 
% and finite extent in the x-direction. Evolution proceeds 
% in discrete time steps, and at each time a check is made 
% to see whether the ball's updated position would be 
% below the floor.  If so, the the y component of the 
% ball's velocity is reversed.
%  
% Finite Difference Approximation: (first order forward 
% differencing,also known as Euler's method).
%  
%   n+1     n 
%  v    =  v  + deltat a
%  
%   n+1     n           n
%  r    =  r  + deltat v
%   
% where the acceleration a is [0, -g] with g constant.
%  
% NOTE: In contrast to the scripts for the pendulum
% simulations, as well as what I will advocate for N-body 
% simulations and others, we do NOT store all the discrete 
% positions and velocities that are generated.  Rather, we 
% simply maintain a 2-element vector for each of the 
% position and velocity, updating the values as we iterate 
% through the time-steps and visually "recording" the 
% values as we go by plotting the ball's position and, 
% optionally, converting all of the "snapshots" to an AVI 
% movie.
% 
% IMPORTANT!!
%
% Please do NOT use the above FDA for term projects 
% involving the solution of dynamical ODEs (e.g. N-body), 
% with the possible exception of initial code development.
% 
%-----------------------------------------------------------
% KEY SIMULATION CONTROL PARAMETERS
%-----------------------------------------------------------
% 
% The following is a list of the key control parameters 
% for the simulation and thus are variables that you can 
% and should experiment with interactively, noting their 
% effect.  It's a good idea to keep track of what you 
% do, analogously to how you are being trained to maintain 
% lab notebooks---except here you're allowed to delete 
% mistakes!).  If you DO maintain a log file, don't call
% it README, since that file should only contain the info
% requested in the handout.
% 
% tmax
% 
%    Final time (length) of simulation.  Depending on the 
%    other parameter settings, and especially when running 
%    the functionized version of the extended version of 
%    this script that you will code, you may need to 
%    adjust tmax to ensure that the simulation time is long 
%    enough for you to study whatever it is you want to 
%    study.  However, when running interactively you can 
%    simply set tmax to a large value, then interrupt the 
%    calculation using, e.g. Ctrl-C, when you've seen enough.
% 
% deltat    
% 
%    Discrete time step (grid spacing, mesh spacing).  This 
%    is THE key numerical control parameter for the 
%    simulation. We expect the simulation to become 
%    increasingly accurate as deltat -> 0, and, as with the 
%    solution of any problem with FDAs we NEED to study the 
%    effect of varying deltat while holding all other 
%    parameters fixed in order to determine how accurate the 
%    simulation is, or more precisely, how accurately we are 
%    solving the equations that define the simulation.
% 
% x0, y0
% 
%    Initial position of the ball. 2 of the 4 parameters that 
%    are REQUIRED in order for the (discrete) equations of 
%    motion to have a solution.  These parameters WILL 
%    influence the overall dynamics, but generally not in a 
%    very interesting way.
% 
% vx0, vy0
% 
%    Initial velocity of the ball.  Specified as a two 
%    element vector in the code.  The other 2 of the 4 
%    parameters needed to yield a solution of the equations 
%    of motion. Changing these will generally more directly 
%    lead to variation in the dynamics, so they are more 
%    interesting to experiment with than positions.
% 
% kappa
% 
%    Coefficient of restitution.  Defines the "bounciness" of 
%    the particle-floor interaction.  Most definitely an 
%    interesting parameter to play with.
%
%-----------------------------------------------------------
% KEY PLOTTING & ANIMATION PARAMETERS
%-----------------------------------------------------------
%
% plotenable
%
%    Set to non-zero/zero to enable/disable plotting
%
%    Disabling plotting will be useful when we use the 
%    functionized version of bounce to perform parameter
%    space surveys ("numerical experiments")
% 
% pausesecs
%
%    The amount of time Matlab will pause and do nothing
%    during each time step of the simulation.  Can be used
%    to control the speed of the interactively displayed
%    animation.  Will have NO EFFECT on the appearance of
%    any AVI movie that is made.
% 
%    Set to 0 for maximum animation speed, something on 
%    the order of a tenth of a second for slow-motion
%    effect.
%
%  avienable
%
%    Set to non-zero/zero to enable/disable recording
%    of an AVI movie.  Enabling recording will slow
%    the real-time animation speed significantly.
%
%    If plotenable is 0, avienable will be set to 0
%    irrespective of what it was on entry to the 
%    function.
%
%  avifilename
%
%    Name for AVI movie file.  Should have an .avi 
%    extension.
%
%  aviframerate
%
%    Presumed/estimated AVI playback speed in frames per
%    second. Used at the start of the simulation to output
%    enough copies of the first frame to ensure that the 
%    ball doesn't move immediately when the movie starts.
%
%-----------------------------------------------------------
% NOTE ABOUT COMMENTING
%
% Clearly, this script is EXTENSIVELY commented, arguably 
% excessively so. You are not expected to comment the 
% scripts and functions that you write for this problem,
% nor those for your term projects, in anything like this 
% sort of detail.  I've been verbose here so that the
% code is self-contained and fully explains both the 
% simulation and the plotting techniques used.
%
% The hope is that all of you will be able to understand 
% it and modify with a minimum of additional help/hints.
%
% NOTE ABOUT VERTICAL WHITESPACE
% 
% Many programming-oriented enterprises, Google for example,
% advocate a programming style that absolutely minimizes
% the number of empty lines.  It is argued that this makes 
% code easier to read and understand, and there can be 
% no dispute that it allows us to see more lines of 
% executable code on the screen. 
% 
% This code obviously violated this dictum, but I believe
% that, especially for programming novices, whitespace improves 
% readability for code that is being perused for the first
% time, particularly if you aren't viewing it with syntax 
% highlighting. 
% 
% However, this is merely  my opinion and you should feel free 
% to comment in any way you see fit.
%-----------------------------------------------------------

%-----------------------------------------------------------
%% SECTION 1. SET UP THE SIMULATION. 
%-----------------------------------------------------------
% Most of this section involves definitions of parameters,
% including quite a few that control aspects of the program
% that don't DIRECTLY relate to the simulation, e.g.
% the color of the ball, the thickness of the line 
% that will represent the floor, the name of the AVI
% file that will be (optionally) created etc.
%-----------------------------------------------------------

%-----------------------------------------------------------
% Parameters defining discrete set of times.
%-----------------------------------------------------------
tmax = 3;
deltat = 0.01;   

%-----------------------------------------------------------
% Coefficient of restitution. Reflected velocity component 
% is a factor of kappa times the incident component.  
%
% Here we set kappa to a value < 1 to compensate for the 
% effective energy gain due to the use of a low accuracy FDA. 
% I.e. if kappa == 1, the ball will bounce higher
% after every reflection (try it!)  Simulate loss of energy
% due to a bounce by setting kappa smaller, or a gain
% ("activating floor) by setting kappa > 1.
%-----------------------------------------------------------
kappa = 0.98;

%-----------------------------------------------------------
% Acceleration of ball.  a is a two-element vector with
% constant x- and y-components. I've tuned g and the initial
% y-velocity to produce a reasonable looking simulation.
%-----------------------------------------------------------
g = -20;
a = [0.0, g];

%-----------------------------------------------------------
% Initial position and velocity of ball. r and v are two
% component vectors.
%-----------------------------------------------------------
x0 = 0.0;
y0 = 1.0;
r = [x0, y0];

vx0 = 0.5;
vy0 = 0.0;
v = [vx0, vy0];

%-----------------------------------------------------------
% Bounce counter.
%-----------------------------------------------------------
nbounce = 0;

%-----------------------------------------------------------
% Set plotenable to non-zero/zero to enable/disable plotting.
%-----------------------------------------------------------
plotenable = 1;

%-----------------------------------------------------------
% Parameter to control speed of animation.  Script execution
% will pause for pausesecs each time a new frame is drawn.
% 
% Setting this parameter to a "largish" value, say 0.1
% (seconds), will produce a slow-motion effect.
% 
% Set it to 0 for maximum animation speed.
%-----------------------------------------------------------
pausesecs = 0.01;

%-----------------------------------------------------------
% Plot attributes defining the appearance  of the ball.
%-----------------------------------------------------------

% Ball has a (marker) size of 15 ...
ballsize = 15;
% ... it's red ...
ballcolor = 'r';
% ... and it's plotted as a circle.
ballmarker = 'o';

%-----------------------------------------------------------
% Define the coordinates for the floor and the properties for 
% the plot command that define the appearance of the floor.  
%
% We store the x and y coordinates of the floor extremities 
% as a 2 x 2 matrix, fxy, with the x coords in the first row 
% and the y coords in the second.  Storing the x and y 
% coords in columns is clearly equivalent. One rationale 
% for being able to access all (both!) x and y values as 
% a row (or column) of a matrix is that it makes the plotting 
% easier (since plot wants VECTORS of x and y coordinates as
% its arguments.  Not a big deal for a geometric object with
% two points, but this approach will clearly generalize to 
% more complex interfaces, for example.
%-----------------------------------------------------------

% Parameters that define the floor position and length.
fxmin = 0.0;
fxmax = 1.0;
fy = 0.0;

% Define the 2 x 2 matrix of floor coordinates in terms of 
% those parameters.
fxy = [[fxmin, fxmax]; [fy, fy]];

%-----------------------------------------------------------
% Plot attributes defining the appearance of the floor.
%-----------------------------------------------------------

% The floor is blue ...
fcolor = 'b';
% ... and we'll draw it using a fairly thick line.
fwidth = 5;

%-----------------------------------------------------------
% Set avienable to non-zero/zero to enable/disable AVI
% video generation.
%
% Note that AVI recording will slow the animation down 
% significantly in general, and if you record many 
% animations in one session, Matlab may start to use a LOT
% of memory, overall performance of the interface may 
% suffer, and you may want to stop and restart Matlab.
%-----------------------------------------------------------
avienable = 1;

% If plotting is disabled, ensure that AVI generation
% is as well
if ~plotenable
   avienable = 0;
end

% Name of avi file.
avifilename = 'bounce_credits.avi';

% Presumed AVI playback rate in frames per second.
aviframerate = 25;

%-----------------------------------------------------------
%% SECTION 2. PERFORM THE SIMULATION.
%-----------------------------------------------------------

%-----------------------------------------------------------
% If AVI creation is enabled, then initialize an avi object 
% (think of it as the AVI file itself).  We absolutely need 
% to do this if we want to record the animation as an 
% AVI file.  The object (handle) that is returned from
% VideoWriter is used subsequently in calls to the functions 
% that record the individual frames, (writeVideo), and 
% close/finalize the AVI file (close).
%-----------------------------------------------------------
if avienable
   aviobj = VideoWriter(avifilename);
   open(aviobj);
end

%-----------------------------------------------------------
% FOR EACH TIME STEP ...
%-----------------------------------------------------------

for t = 0 : deltat : tmax
   % BEGIN TIME STEP
    
   if plotenable
      %%---------------------
      %% BEGIN Graphics section
      %%---------------------

      % The simple animation technique used here dictates that 
      % at every time step we erase the plotting window and 
      % redraw EVERYTHING, including the axes, labels, title, 
      % floor AND, of course, the ball.

      % Clear figure
      clf;

      % Don't erase figure after each plot command.
      hold on;

      % Define plotting area, using a 1:1 aspect ratio for the 
      % plotted region ...
      axis square;
      % ... boxed axes ...
      box on;
      % ... and including a 15%-width "border" around the unit square.
      dlim = 0.15;
      xlim([-dlim, 1 + dlim]);
      ylim([-dlim, 1 + dlim]);

      % Make and display title.  Note the use of sprintf to create
      % a string containing formatted values of the time step,
      % the time, and the number of bounces.
      titlestr = sprintf('Step: %d   Time: %.1f  Number of bounces: %d', ...
         fix(t / deltat), t, nbounce);

      % This illustrates how to set the font size/weight and 
      % RGB color of text. The RGB format, a three element row vector 
      % of R, G and % B values between 0.0 and 1.0 inclusive, can 
      % also be used with other plot elements, such as lines and markers.
      title(titlestr, 'FontSize', 16, 'FontWeight', 'bold', ...
         'Color', [0.25, 0.42, 0.31]);

      % Draw the floor.
      plot(fxy(1,:), fxy(2,:), 'Color', fcolor, 'LineWidth', fwidth);

      % Draw the ball. Set MarkerEdgeColor and MarkerFaceColor to
      % the same color to produce a filled circle.
      plot(r(1), r(2), 'Marker', ballmarker, 'MarkerSize', ballsize, ...
         'MarkerEdgeColor', ballcolor, 'MarkerFaceColor', ballcolor);

      % Force update of figure window.  This statement can be important 
      % when using the simple animation technique demonstrated here: it 
      % forces any pending plotting operations that might have been 
      % "buffered" to be executed, i.e. it ensures that now that we 
      % have plotted everything we need to plot at any given time step,
      % we will actually see what we have drawn.
      drawnow;

      % Record video frame if AVI recording is enabled.
      if avienable
         % Is this the first time step?
         if t == 0
            % Define framecount so that multiple copies of the first 
            % frame, corresponding to about 5 seconds of the video,
            % are recorded. If the movie comes up running, the 
            % viewer will then have a bit of time to absorb the 
            % physical setup.
            framecount = 5 * aviframerate ;
         else
            % For all other times, record a single frame. 
            framecount = 1;
         end
         for iframe = 1 : framecount
            writeVideo(aviobj, getframe(gcf));
         end
      end

      % Pause execution.
      pause(pausesecs);
      %%---------------------
      %% END Graphics section
      %%---------------------
   end

   %%--------------------------
   %% Particle dynamics section
   %%--------------------------

   % Note how few lines of code are needed to do the actual 
   % "simulating" --- about 10 out of the 75 or so non-comment,
   % non-whitespace lines.  This is typical of a scientific
   % computing program.  The "infrastructure" that we need to 
   % control and analyze the calculation often completely 
   % dominates (in terms of length of code) the part of the 
   % program that does most of the calculational work.

   % Note that r, v, dr and a are all two-component row vectors
   % (x and y components), and as noted in the header comments,
   % that we do NOT store all of the positions and velocities
   % as they are computed.

   % Update the velocity.
   v = v + deltat * a;

   % Compute provisional change in ball's position.
   dr = deltat * v;

   % Check for an intersection with the floor. There are 
   % two conditions: The first comparison determines whether the 
   % particle's provisional y-coordinate is  below the level 
   % of the floor, while the second and third comparisons determine 
   % whether the particle's previous horizontal location positions 
   % it over the floor, i.e so that it hasn't bounced past the 
   % right edge of the floor.  
   % 
   % Recall that the value of the first index into the 
   % 2D array (matrix) fxy, tells us whether we are referring 
   % to an x (1) or y (2) coordinate,while the second tells us 
   % whether we are referring to the left endpoint (1) or 
   % right (2).  Also, since the floor is horizontal, we have 
   % fxy(2,1) = fxy(2,2), so it doesn't matter which we use in 
   % the first comparison.

   if ( r(2) + dr(2) < fxy(2,1) )  &  ( fxy(1,1) <= r(1) & r(1) <= fxy(1,2) )

      %-------------------------------------------------------
      % Provisional particle position IS below the floor.
      %-------------------------------------------------------

      % Reverse sign of v_y, and scale it by the coefficient of 
      % restitution.
      v(2) = - kappa * v(2);

      % Don't change the particle's y coordinate.
      dr(2) = 0;

      % Update the bounce count.
      nbounce = nbounce + 1;

   end 

   % Update the position.
   r = r + dr;

   % END TIME STEP
end

% END SIMULATION

% If we're making a video, record the credit frames and finalize 
% the recording.

if avienable
   % Record 3 seconds of final frame.
   for iframe = 1 : 3 * aviframerate 
      writeVideo(aviobj, getframe(gcf));
   end
   text(0.5, 0.65, 'Filmed on location in Hennings 205.', ...
      'FontSize', 16, 'FontWeight', 'bold', 'Color', [1, 0, 0], ...
      'HorizontalAlignment', 'center');
   % Record 5 seconds of first credit frame.
   for iframe = 1 : 5 * aviframerate 
      writeVideo(aviobj, getframe(gcf));
   end
   text(0.5, 0.35, ...
      'No balls were harmed in the making of this video.', ...
      'FontSize', 12, 'FontWeight', 'bold', 'Color', [1, 0, 0], ...
      'HorizontalAlignment', 'center');
   % Record 5 seconds of second credit frame.
   for iframe = 1 : 5 * aviframerate
      writeVideo(aviobj, getframe(gcf));
   end

   % Finalize recording and write a diagnostic message that 
   % a movie file was created.
   close(aviobj);
   fprintf('Created video file: %s\n', avifilename);

end
% END OF SCRIPT
