%newbounce - a new attempt to model the 'bouncing ball' problem'
%It still has some kinks, but is FAR more well behaved than the 
%previous version - I've moved the contact detection out of the
%original loop. I think this made it easier to see what was
%happening - and finally got the coefficient of restitution to
%take energy OUT of the system.

%Initialisation - quick version - limiting some variables for test speed
e = input('Enter e ::  ');
y_init = input('Starting height? ::  ');
k = input('Spring constant (recommend about 100)? ::  ');
step = 1/2000;
N = 30;
r = 0.2;
deltamax = r/100;
g = 9.81;
m = 1;

y = y_init;
y_dot =0;

delta = y-r;

gravity_force = (-m*g);
contact_force = 0;	%assuming we start off the ground
contact = 0 		%contact detection variable/trigger

Results = zeros(1,(N*1/step));
Delta = zeros(1,(N*1/step));
Y_dot = zeros(1,(N*1/step));
Contact = zeros(1,(N*1/step));
Contact_force = zeros(1,(N*1/step));

disp('Press a key to begin');

for n = 0:step:N

	%Contact detection
	if y <= r
		contact = 1;
	elseif y > r 
		contact = 0;
	end

	%No contact - falling under gravity	
	if contact == 0
		y_dot = y_dot + (gravity_force * (step/m));
		y = y + (step * y_dot);
		contact_force = 0;
		Results(round(n/step + 1)) = y-r;
		
	%Contact - impact and restitution
	elseif contact == 1
		disp('Boing');
			
			if y_dot <= 0					%This should produce
				k_live = k;					%a kind of hysteresis
			elseif y_dot > 0				%using the coefficient
				k_live = (k * (e^2));		%of resititution
			end
			
		
		contact_force = k_live * delta;
		y_dot = y_dot + (contact_force * (step/m));
		y = y + (y_dot * step);
		Results(round(n/step + 1)) = y-r;
	end
	Delta(round(n/step + 1)) = delta;
	Y_dot(round(n/step + 1)) = y_dot;
	Contact(round(n/step + 1)) = contact;
	Contact_force(round(n/step + 1)) = contact_force;
end
hold on
plot(Results);