function [q_n,p1,p0]=hybrid_mc(q_0,hmc_parms,W_d)

parms     =  hmc_parms{1};
bbk_parms =  hmc_parms{2};
nsteps    =  hmc_parms{3};

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

q_orig = q_0;
%select momentum at random
sig    = sqrt(diag(M));
v_0    = randn(2,1).*sig;
KE_0   = (M*v_0)'*v_0/2;
U_0    = -log(e_V( q_0(1),q_0(2),parms ) );

for i=1:nsteps
   [q_n,v_n]=BBK_r2(q_0,v_0,bbk_parms,parms);
   q_0 = q_n ;  v_0 = v_n ;
   W_d=
end

KE = (M*v_0)'*v_0/2;
U  = -log(e_V(q_0(1),q_0(2),parms));

p0=exp(-KE_0 - U_0);
p1=exp(-KE   - U  );
% ratio=p1/p0;
% 
% if (ratio>=1)
%     accept=1;
% else
%     if (rand<ratio)
%         accept=1;
%     else
%         accept=0;
%     end
% end
% 
% if (accept==0)
%     q=q_orig;
% else
%     q=q_0;
% end