function [pos_final,vel_final] = BR_bounded_motion3D(box,p0,v0,del_t,...
    faces,normals)
%

V_VEL = [-10,0];
DAMPING = 0.9;
MAX_VEL = 2000;

x_min = min(box(:,1));
x_max = max(box(:,1));
y_min = min(box(:,2));
y_max = max(box(:,2));
z_min = min(box(:,3));
z_max = max(box(:,3));

pos = p0;
vel = v0;
speed = norm(vel);
p_new = pos + del_t*vel;
dir = vel/norm(vel);
dist2go = norm(p_new-pos);
while dist2go>0
    if p_new(1)>=x_min&p_new(1)<=x_max...
            &p_new(2)>=y_min&p_new(2)<=y_max...
            &p_new(3)>=z_min&p_new(3)<=z_max
        pos_final = p_new;
        vel_final = vel;
        return
    end
    [ip,dist2int,normal] = BR_int_box(faces,normals,pos,p_new);
    pos = ip;
    v = cross(-dir,normal);
    v = v/norm(v);
    theta = acos(dot(-dir,normal));
    R = BR_gen_R(v,2*theta);
    dir = (R*dir')';
    if dot(dir,normal)<0
        dir = -dir;
    end
    vel = dir*speed;
    dist2go = max(0,dist2go-dist2int);
    p_new = pos + vel*del_t;
end
pos_final = pos;
vel_final = dir;
if norm(vel)>MAX_VEL
    vel = MAX_VEL*vel/norm(vel);
end
