%flash_def.m
%   pot1.m
%   pot2.m
%   fk12.m
%   fk21.m
% This function computes the mean velocity and effective diffusion 
%coeficient of a particle moving in a flashing potential. 

function [vel,Deff] = flash_deff(ngrid, lp, D, kT, ppot1,...
    ppot2, pk12, pk21, pot1, pot2, fk12, fk21)


dx = lp/ngrid; %grid spacing
%set up the grid
xpos = dx/2 +[0:1:ngrid-1]*dx;

gam = D/dx^2; %Diffusive transition rate

%evaluate chemical transition rates 
k12 = feval(fk12, xpos, pk12);
k21 = feval(fk21, xpos, pk21);
 
%Compute the delta phi's
dphi1 = feval(pot1,[xpos(2:ngrid),xpos(ngrid)+dx],ppot1) ... 
       - feval(pot1,xpos,ppot1);
dphi2 = feval(pot2, [xpos(2:ngrid),xpos(ngrid)+dx],ppot2) ... 
       - feval(pot2,xpos,ppot2);
   
%compute the spatial transistion rates
for i=1:ngrid
    s = dphi1(i)/kT;
  if abs(s) > 1.5e-3
   F1(i) = gam*s/(exp(s) - 1);
   B1(i) = gam*(-s)/(exp(-s) - 1);
  else
   F1(i) = gam/(s*(s*(s/24+1/6)+1/2)+1);
   B1(i) = gam/(1-s*(1/2-s*(1/6-s)));
  end
   s = dphi2(i)/kT;
  if abs(s) > 1e-6
   F2(i) = gam*s/(exp(s) - 1);
   B2(i) = gam*(-s)/(exp(-s) - 1);
  else
   F2(i) = gam/(s*(s*(s/24+1/6)+1/2)+1);
   B2(i) = gam/(1-s*(1/2-s*(1/6-s)));
  end
end

%set up the matrices
L1=diag(-(F1+[B1(ngrid),B1(1:ngrid-1)])) ...
   +diag(F1(1:ngrid-1),-1)+diag(B1(1:ngrid-1),1);
L1_p = zeros(ngrid, ngrid);
L1_p(1,ngrid) = F1(ngrid);
L1_m = zeros(ngrid, ngrid);
L1_m(ngrid,1) = B1(ngrid);

L2=diag(-(F2+[B2(ngrid),B2(1:ngrid-1)])) ...
   +diag(F2(1:ngrid-1),-1)+diag(B2(1:ngrid-1),1);
L2_p = zeros(ngrid, ngrid);
L2_p(1,ngrid) = F2(ngrid);
L2_m = zeros(ngrid, ngrid);
L2_m(ngrid,1) = B2(ngrid);

M1 = L1+L1_m+L1_p;
M2 = L2+L2_m+L2_p;

M = zeros(2*ngrid,2*ngrid);
L_p = zeros(2*ngrid,2*ngrid);
L_m = zeros(2*ngrid,2*ngrid);

M(1:ngrid,1:ngrid) = M1;
M(ngrid+1:2*ngrid,ngrid+1:2*ngrid) = M2;
M = M-diag([k12,k21])...
    +diag(k21,ngrid)+diag(k12,-ngrid);
L_p(1:ngrid,1:ngrid) = L1_p;
L_p(ngrid+1:2*ngrid,ngrid+1:2*ngrid) = L2_p;
L_m(1:ngrid,1:ngrid) = L1_m;
L_m(ngrid+1:2*ngrid,ngrid+1:2*ngrid) = L2_m;



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

%the left zero eigenvector
l0 = ones(1,2*ngrid);

%lamp = l0*[L_p-L_m]*ps;
vel = l0*[L_p-L_m]*ps;

rhs =  (l0*(L_p-L_m)*ps)*ps - (L_p-L_m)*ps;
rhs(2*ngrid,1) = 0;

rp = Mp\rhs;

Deff =  lp^2/2 *(l0*(L_p+L_m)*ps +2*l0*(L_p-L_m)*rp);

%lampp = l0*(L_p+L_m)*ps +2*l0*(L_p-L_m)*rp;


%vel = lp*lamp;
%Deff = lp^2/2 *(lampp);