%% Two-step system simulation - AJP
% Marie Lopez del Puerto
 
% 1. Inputs

% Width of box:
W = 0.2; % meters
% Length of lower level 
L0 = 0.2; % meters
% Length of upper level:
L1 = 0.4; % meters
% Height difference between levels:
h = 0.01;
% Gravitational acceleration:
g = 9.81; % m/s^2
% Time-step:
dt = 0.01; % seconds
% Number of time-steps:
N = 1e5; 
% Maximum velocity of 'kicks':
vmax = 1;

%%

t0 = 0; % counter

% 2. Generate a random initial position vector:

x = [(L0+L1)*rand W*rand];  % rand generates a random number between 0 and 1.

% 3. Generate a random initial velocity vector:

r1=rand;
r2=rand;
v = [r2*cos(2*pi*r1) r2*sin(2*pi*r1)]; % in m/s.

% 4. Loop:

for i = 2:N+1
    
    xp = x; % save previous value of x
    
    % a) Modify position:
    
    x = x + v*dt; 
    
    % b) If statements to check for wollisions with walls:
    
    KE = dot(v,v)/2; % Kinetic energy divided by mass
    PE = g*h; % Potential energy divided by mass
    
    if x(1) < 0 % left wall
        v(1) = -v(1);
        x(1) = -x(1);
    end
    
    if x(1) > L0 && xp(1) < L0  % step going up
        if KE < PE
            v(1) = -v(1);
            x(1) = 2*L0-x(1);
        else
            ratio = sqrt(1-PE/KE); % from energy conservation
            v = ratio*v;
        end
    end
    
    if x(1) < L0 && xp(1) > L0  % step going down
            ratio = sqrt(1+PE/KE); % from energy conservation
            v = v./ratio;
    end
    
    if x(1) > L0 + L1 % right wall
        v(1) = -v(1);
        x(1) = 2*(L0+L1)-x(1);
    end
    
    if x(2) < 0 % bottom wall
        v(2) = -v(2);
        x(2) = -x(2);
    elseif x(2) > L0 % top wall
        v(2) = -v(2);
        x(2) = 2*L0-x(2);
    end

    if x(1) < L0
        t0 = t0 + 1;
    end
    
    % e) Every X time-steps, randomy modify v:
        
    if mod(i,100) == 0                     
        r1=rand;                            
        r2=rand;           
        v = vmax*[r2*cos(2*pi*r1) r2*sin(2*pi*r1)]; 
    end
    
end

t1 = N-t0;


%% Display results:

disp(['Ratio of areas, a1/a0 = ', num2str(L1/L0)]);
disp(['Step height, h = ', num2str(h)])
disp(['Number of time steps spent in lower level, t0 = ', num2str(t0)]);
disp(['Number of time steps spent in upper level, t1 = ', num2str(t1)]);
disp(['t1/t0 = ', num2str(t1/t0)]);