function [q,v]=BBK(q_in,v_in,bbk_parms,parms);

M=bbk_parms{1};
dt=bbk_parms{2};
gamma=bbk_parms{3};

%sig=sqrt( 2*gamma*diag(M) );
%keyboard;

%1st start with Verlet:

dp=delta_p(q_in,v_in,gamma,bbk_parms,parms);
v_mid = (1-0.5*gamma*dt)*v_in + M^-1*dt/2*dp;
q     = q_in + v_mid*dt;

dp=delta_p(q,v_mid,gamma,bbk_parms,parms);

v=v_mid + M^-1*dt/2*dp;
v=v/( (1 + 0.5*gamma*dt) );
end

function dp=delta_p(q_in,v_in,gamma,bbk_parms,parms)
M=bbk_parms{1};
dt=bbk_parms{2};
gamma=bbk_parms{3};
sig=sqrt( 2*gamma*diag(M)/dt );
GV=grad_V(q_in(1),q_in(2),parms);
if(gamma>0) 
    R=randn(2,1).*sig;
    dp = -GV + R;  
else
    dp = -GV;
end
end
