%%%
% Dear students:
% This is a little program to demostrate the motion of a particle subject to gravitational acceleration
% The initial conditions of the particle are given by its initial velocity v0 (in m/s) and 
% by the initial angle theta(in degree).
%
% Although you use MathLab in the framework of our STEPS program in math and engineering, you probably
% are not familiar with this kind of program. However, I encourage you to play a little bit around. You
% will also understand better the meaning of integrating a differential equation. What happens, if you
% use bigger time steps, i.e. a bigger Deltat ?
%
%%%
function  vdpode(MU)

x0 = 0;             % initial position, x-component (in m)
y0 = 0;             % initial position, y-component (in m)
v0 = 10;            % initial velocity of the particle (in m/s)
theta = 40;         % intial angle of the particle (in degree)
ax = 0;             % accerleration in x direction;   
ay = -9.8;          % accerleration in y direction (in m/s^2);           
x = zeros(1000,1);     % x(t)
y = zeros(1000,1);     % y(t)
vx = zeros(1000,1);    % vx(t)
vy = zeros(1000,1);    % vy(t)
t= zeros(1000,1);      % t(i) = i*Deltat
Deltat = 0.005;        % integration step of time (in second) 

i = 1;
vx(i) = v0*cos(theta/180*pi);
vy(i) = v0*sin(theta/180*pi);
x(i) = x0;
y(i) = y0;
t(i) = 0;
while y(i) >= -0.0000
    i = i + 1;
    t(i) = t(i-1) + Deltat;
    x(i) = x(i-1) + vx(i-1) * Deltat;
    y(i) = y(i-1) + vy(i-1) * Deltat;
    vx(i) = vx(i-1) + ax * Deltat;
    vy(i) = vy(i-1) + ay * Deltat;
    figure(1);
    plot(t(1:i),x(1:i),'ok-');
    title(['x position of particle']);
    xlabel('time t / second');
    ylabel('X(t) / meter');
    xlim([0 1.5]);
    ylim([0 12]);
    figure(2);
    plot(t(1:i),y(1:i),'or-');
    title(['y position of particle']);
    xlabel('time t / second');
    ylabel('Y(t) / meter');
    xlim([0 1.5]);
    ylim([0 2.5]);
    figure(3);
    plot(x(1:i),y(1:i),'or-');
    title(['x-y position of particle']);
    xlabel('x(t) / meter');
    ylabel('Y(t) / meter');
    xlim([0 12]);
    ylim([0 2.5]);
    pause(0.01);
end