% avgv_deff.m
%    pot1.m 
%This program computes the mean velocity and effective diffusion 
%coeficient of a particle moving 1-D tilted-periodic potential. 
%The user must supply the function pot1.m that evaluates the 
%potential.

%Function inputs
%   Name            Description
%   ngrid           the number of grid points
%   lp -            the length of one period of the potential
%   D               the diffusion coefficient of the particle
%   kT              the Boltzmann constant times the absolute temperature
%   ppot1          a vector containing the parameters needed by the
%                   function pot1.
%   pot1(x,ppot1)   user supplied function that evaluates the potential
%                   at postion x using the parameters in the vector ppot1.

function [v,d] = avgv_deff(ngrid, lp, D, kT, ppot1, pot1)


dx = lp/ngrid;  %grid spacing
xpos = dx/2 +[0:1:ngrid-1]*dx;  %set up the grid
gam = D/dx^2;  %diffusive transition rate

%compute Delta phi
dphi = feval(pot1,[xpos(2:ngrid),xpos(ngrid)+dx],ppot1) ... 
     - feval(pot1,xpos,ppot1); 
   
%set up the transition rates
for i=1:ngrid
    s = dphi(i)/kT;
  if abs(s) > 1.5e-3
   F(i) = gam*s/(exp(s) - 1);
   B(i) = gam*(-s)/(exp(-s) - 1);
  else
   F(i) = gam/(s*(s*(s/24+1/6)+1/2)+1);
   B(i) = gam/(1-s*(1/2-s*(1/6-s)));
  end
end

%Set up the matrices
L=diag(-(F+[B(ngrid),B(1:ngrid-1)]))+diag(F(1:ngrid-1),-1)+diag(B(1:ngrid-1),1);
L_p = zeros(ngrid, ngrid);
L_p(1,ngrid) = F(ngrid);
L_m = zeros(ngrid, ngrid);
L_m(ngrid,1) = B(ngrid);
M = L+L_m+L_p;

%Solve for the steady-state probabilities
Mp = M;
Mp(ngrid,:) = ones(1,ngrid);
f = zeros(ngrid,1);
f(ngrid,1) = 1;
ps = Mp\f;

%The left zero eigenvector
l0 = ones(1,ngrid);

%compute the velocity
v = lp*l0*[L_p-L_m]*ps;

%compute the effective diffusion coefficient
rhs =  (l0*(L_p-L_m)*ps)*ps - (L_p-L_m)*ps;
rhs(ngrid,1) = 0;
rp = Mp\rhs;
d = lp^2/2*(l0*(L_p+L_m)*ps +2*l0*(L_p-L_m)*rp);