%samp_paths.m
%   pot1.m
%   pot2.m
%   fk12.m
%   fk21.m
% This function generates sample paths of the flashing ratchet.


function [tvals,xvals,svals] = samp_paths(xin,sin, tfinal, nsamp, dx, D, kT, ppot1,...
    ppot2, pk12, pk21, pot1, pot2, fk12, fk21)

gam = D/dx^2; 

time = 0;
x = xin;
state = sin;

%sample period
stime = tfinal/nsamp;

dtime = stime;

%ouput vectors
tvals = zeros(1,nsamp+1);
xvals = zeros(1,nsamp+1);
svals = zeros(1,nsamp+1);

tvals(1) = time;
xvals(1) = x;
sval(1) = state;

i = 1;

while time < tfinal
     
    if state == 1
        k = feval(fk12,x,pk12);
         dphif = feval(pot1,x+dx,ppot1) ... 
              - feval(pot1,x,ppot1);
         dphib = feval(pot1,x,ppot1) ... 
           - feval(pot1,x-dx,ppot1);
    else
        k = feval(fk21, x,pk21);
        dphif = feval(pot2,x+dx,ppot2) ... 
              - feval(pot2,x,ppot2);
        dphib = feval(pot2, x,ppot2) ... 
             - feval(pot2,x-dx,ppot2);
    end
 
    s = dphif/kT;
    if abs(s) > 1.5e-3
        F = gam*s/(exp(s) - 1);
     else
         F = gam/(s*(s*(s/24+1/6)+1/2)+1);
    end
  
    s = dphib/kT;
    if abs(s) > 1.5e-3
        B = gam*(-s)/(exp(-s) - 1);
    else
        B = gam/(1-s*(1/2-s*(1/6-s)));
    end
  
    rate = k+F+B;
    dt = -1/rate*log(rand);
    time = time+dt; 
    p = rand;
    
    if p<k/rate
        state = abs(state-1);
    elseif k/rate<= p & p < (k+F)/rate
        x = x+dx;
    else
        x = x-dx;
    end
    
    if time > dtime
        dtime = dtime+stime;
        svals(i+1) = state;
        xvals(i+1) = x;
        tvals(i+1) = time;
        i = i+1;
    end
end