Bouncing ball's height and velocity!

I want to make a graph when the ball drop from the height(20m). I insert what I wanted! I want to know what I do wrong.

clear all
h0 = 20;      
v = 0;      
g = 10;         
t=0;
dt = 0.01;     
rho = 0.75;     
tau = 0.10 ;    
hmax = h0   ;   
h = h0;
hstop = 0.01;   
freefall = 1; 
t_last = -sqrt(2*h0/g); 
vmax = sqrt(2 * hmax * g);
H = [];
T = [];
while(hmax > hstop)
    if(freefall==1)
        hnew = h + v*dt - 0.5*g*dt*dt;
    if(hnew<0)
         t = t_last + 2*vmax;
         freefall = 0;
         t_last = t + tau;
         h = 0;
    else
      t = t + dt;
      v = v - g*dt;
      h = hnew;
    end
    else
      t = t + tau;
      vmax = vmax * rho;
      v = vmax;
      freefall = 1;
      h = 0;
    end
    hmax = 0.5*vmax*vmax/g;
  H.append(h);
  T.append(t);
end 
  
%% Simulation
plot(time, height, 'r.', 'MarkerSize', 50);
axis([-2, 20, 0 25]); grid;
xlabel('ball position X [m]')
ylabel('ball position Y [m]')
title('TaengTaeng Ball')
drawnow

Answer 1

[Image Analyst]

Well there was a lot that was wrong with that code. Here is the fix:

% Demo to simulate a bouncing ball.
clc;    % Clear the command window.
fprintf('Beginning to run %s.m.\n', mfilename);
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

h0 = 20;	% Initial drop height in meters.
v = 0;		% Initial y velocity in m/sec.
g = 9.8;	% Gravitational acceleration in m/s^2;
t = 0;		% Initial time when dropped
dt = 0.01;	% Delta time in seconds.
rho = 0.75;	% Velocity reduction factor.  Velocity reduces this much after a bounce.
peakHeight = h0;	% Initial drop height in meters.
h = h0;		% Instantaneous height.
hstop = 0.01;	% Height at which if the peak height after a bounce is less than this, stop the simulation.
% Preallocate arrays for time and height.  Make them plenty large - we will crop to the final size later.
T = 0 : dt : 1000;
H = zeros(1, length(T));
% Setup a failsafe.  Don't do more than this number of iterations or else we might have an infinite loop.  This will prevent that.
maxIterations = 100000;		
loopCounter = 1;
while (peakHeight > hstop) && (loopCounter < maxIterations)
	% Compute new height.
	hNew = h + v * dt - 0.5 * g * dt ^ 2;
	% 	fprintf('After iteration %d, time %f, hmax = %f.\n', loopCounter, T(loopCounter), hNew);
	if(hNew<0)
		% Ball hit the ground.
		
		% Find index of last time h was 0
		lastBounceIndex = find(H(1 : loopCounter-1) == 0, 1, 'last');
		if isempty(lastBounceIndex)
			% If it hasn't bounced yet, start looking from the beginning.
			lastBounceIndex = 1;
		end
		% Compute the greatest height since the last bounce, or the initial release.
		[peakHeight, index] = max(H(lastBounceIndex : end)); % Record height
		% Find time when it was at that height.
		tMax = T(index + lastBounceIndex - 1);
		plot(tMax, peakHeight, 'b+', 'MarkerSize', 18, 'LineWidth', 2);
		hold on;
		fprintf('After iteration %d, time %f, hmax = %f.\n', loopCounter, tMax, peakHeight);
		
		% Reflect it up.  For example, if at this time,
		% the ball was going to be at -4 (with no ground in the way)
		% Now, after bouncing, it would be at +4 above the ground.
		h = 0; %abs(hNew);
		v = -v * rho;
		
	else
		% Ball is falling or rising.
		v = v - g*dt;
		h = hNew;
	end
	H(loopCounter) = h;
	loopCounter = loopCounter + 1;
end
% Crop times.
T = T(1 : loopCounter - 1);
H = H(1 : loopCounter - 1);
% Plot the trajectory.
plot(T, H, 'r.', 'MarkerSize', 5);
grid on;
xlabel('Time [seconds]', 'FontSize', fontSize)
ylabel('Ball Position Y [m]', 'FontSize', fontSize)
title('Bouncing Ball', 'FontSize', fontSize)
% Maximize window
g = gcf;
g.WindowState = 'maximized';
fprintf('Done running %s.m.\n', mfilename);

[희현 김]

Can you change time-axis to Ball Position X when Ball has a 0.5m/s velocity??

[Image Analyst]

Sure. Did you try it? It's a really easy change to make:

% Demo to simulate a bouncing ball.
clc;    % Clear the command window.
fprintf('Beginning to run %s.m.\n', mfilename);
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

h0 = 20;	% Initial drop height in meters.
vy = 0;		% Initial y velocity in m/sec.
vx = 0.5;		% Initial x velocity in m/sec.
g = 9.8;	% Gravitational acceleration in m/s^2;
t = 0;		% Initial time when dropped
dt = 0.01;	% Delta time in seconds.
rho = 0.75;	% Velocity reduction factor.  Velocity reduces this much after a bounce.
peakHeight = h0;	% Initial drop height in meters.
h = h0;		% Instantaneous height.
hstop = 0.01;	% Height at which if the peak height after a bounce is less than this, stop the simulation.
% Preallocate arrays for time and height.  Make them plenty large - we will crop to the final size later.
T = 0 : dt : 1000;
H = zeros(1, length(T));
x = vx * T;

% Setup a failsafe.  Don't do more than this number of iterations or else we might have an infinite loop.  This will prevent that.
maxIterations = 100000;		
loopCounter = 1;
while (peakHeight > hstop) && (loopCounter < maxIterations)
	% Compute new height.
	hNew = h + vy * dt - 0.5 * g * dt ^ 2;
	% 	fprintf('After iteration %d, time %f, hmax = %f.\n', loopCounter, T(loopCounter), hNew);
	if(hNew<0)
		% Ball hit the ground.
		
		% Find index of last time h was 0
		lastBounceIndex = find(H(1 : loopCounter-1) == 0, 1, 'last');
		if isempty(lastBounceIndex)
			% If it hasn't bounced yet, start looking from the beginning.
			lastBounceIndex = 1;
		end
		% Compute the greatest height since the last bounce, or the initial release.
		[peakHeight, index] = max(H(lastBounceIndex : end)); % Record height
		% Find time when it was at that height.
		tMax = T(index + lastBounceIndex - 1);
		xMax = x(index + lastBounceIndex - 1);
		plot(xMax, peakHeight, 'b+', 'MarkerSize', 18, 'LineWidth', 2);
		hold on;
		fprintf('After iteration %d, time %f, hmax = %f, xMax = %f.\n', loopCounter, tMax, peakHeight, xMax);
		
		% Reflect it up.  For example, if at this time,
		% the ball was going to be at -4 (with no ground in the way)
		% Now, after bouncing, it would be at +4 above the ground.
		h = 0; %abs(hNew);
		vy = -vy * rho;
		
	else
		% Ball is falling or rising.
		vy = vy - g*dt;
		h = hNew;
	end
	H(loopCounter) = h;
	loopCounter = loopCounter + 1;
end
% Crop times.
T = T(1 : loopCounter - 1);
x = vx * T;
H = H(1 : loopCounter - 1);
% Plot the trajectory.
plot(x, H, 'r.', 'MarkerSize', 5);
grid on;
xlabel('X [meters]', 'FontSize', fontSize)
ylabel('Ball Position Y [m]', 'FontSize', fontSize)
title('Bouncing Ball', 'FontSize', fontSize)
% Maximize window
g = gcf;
g.WindowState = 'maximized';
fprintf('Done running %s.m.\n', mfilename);

[희현 김]

Your answer is so helpful!!!! I can do what I want. Thanks!!!

[dan p]

Is it possible to rewrite the code for "for loop" istead of "while"?

[Image Analyst]

@dan p, of course you can convert the while loop into a for loop. You can either iterate a fixed number of bounces, or have it be like a million bounces and just break whenever you want, like when the height is less than some small value. I'm sure you can do it if you just try.

[dan p]

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.

h0 = 1;	            % Initial drop height in meters.
v = 0;		        % Initial y velocity in m/sec.
g = 9.8;	        % Gravitational acceleration in m/s^2;
t = 0;		        % Initial time when dropped
dt = 0.01;	        % Delta time in seconds.
rho = 0.8;	        % Velocity reduction factor.  Velocity reduces this much after a bounce.
peakHeight = h0;	% Initial drop height in meters.
h = h0;		        % Instantaneous height.
hstop = 0.03;	% Height at which if the peak height after a bounce is less than this, stop the simulation.
ntimes = 10000

ntimes = 10000

% Preallocate arrays for time and height.  Make them plenty large - we will crop to the final size later.
T = 0 : dt : 1000;
H = zeros(1, length(T));
	
% loopCounter 
i = 1;
for (i = 1:ntimes)
	% Compute new height.
	hNew = h + v * dt - 0.5 * g * dt ^ 2;
    
	if(hNew<0)
		% Ball hit the ground.
		
		% Find index of last time h was 0
		lastBounceIndex = find(H(1 : i-1) == 0, 1, 'last');
		if isempty(lastBounceIndex)
			% If it hasn't bounced yet, start looking from the beginning.
			lastBounceIndex = 1;
        end
		
		h = 0; 
		v = -v * rho;
		
	else
		% Ball is falling or rising.
		v = v - g*dt;
		h = hNew;
	end
	H(i) = h;
	i = i + 1;
end

% Crop times.
T = T(1 : i - 1);
H = H(1 : i - 1);

% Plot the trajectory.
plot(T, H);
axis([0, 3, 0 1]);
grid on;
xlabel('Time [seconds]')
ylabel('Ball Position [m]')
title('Bouncing Ball')

My code is for 1m but it worked with "for loop" big thanks :)))

[Image Analyst]

@dan p, at the bottom of the for loop, if you wanted to stop when h < hstop, you'd put

if c <= hstop
    break;
end

It won't get to a million iterations then.

[dan p]

@Image Analyst Maybe stupid question, but what is c?

[Image Analyst]

@dan p, sorry, it should have been the current height, h

if h <= hstop
    break;
end

By the way it's not recommended to use i as a variable since it's also used as the imaginary constant sqrt(-1).

Answer 2

[Image Analyst]

You never defined the "time" vector. Also, you can't do this:

  H.append(h);
  T.append(t);

since H and T are null and don't have an append method.

Also there is a glaring lack of comments.

Did you try to translate this from another language, or did you write this from scratch in MATLAB?