x0=-5;
y0=0;
z0=200;
v0x=2;
v0y=4;
v0z=0;

m=1;
k=0.8;
g=10;

p0=[x0, y0, z0];

points=[];
points=[points ; p0];

pointsToDraw=[];
pointsToDraw=[pointsToDraw ; p0];

tS=0.001;
t =0;

for i=1:1:5
    a=-v0x^2 - v0y^2 - 0.5 * g;
    b=-2*x0*v0x - 2*y0*v0y + v0z;
    c=-x0^2 - y0^2 + z0;

    delta=b^2-4*a*c;

    if(delta>=0)
        t1=(-b - sqrt(delta))/(2*a);
        t2=(-b + sqrt(delta))/(2*a);

        t=max(t1,t2);

        x1=x0+v0x*t;
        y1=y0+v0y+t;
        z1=z0+v0z+t - 0.5*g*t^2;
        points = [points; x1 y1 z1];
        pointsToDraw=[pointsToDraw; x1 y1 z1];

        x0=x1;
        y0=y1;
        z0=z1;


        N = [-2*x1,-2*y1, 1];

        n = N/norm(N);
        vb = [v0x v0y v0z-g*t];
        va = sqrt(k) * (vb - (2 * (dot(vb, n) * n)));

        v0x = va(1);
        v0y = va(2);
        v0z = va(3);

        display(t);
    end
end

% points;
% velocities;
% [X,Y] = meshgrid(-20:1:25);