% flash_int.m
%   pot1.m
%   pot2.m
%   fk12.m
%   fk21.m
%This function solves the diffusion equation for a two-state
%flashing ratchet.

function [rho,xpos] = flash_int(p,ngrid, lp, D, kT, tfinal, dt, ppot1,...
    ppot2, pk12, pk21, pot1, pot2, fk12, fk21)


nsteps = ceil(tfinal/dt); %number of steps to iterate
dx = lp/ngrid; %grid size

%Evaluate the grid
xpos = dx/2 +[0:1:ngrid-1]*dx;

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

%Evaluate the 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 overall transition matrix
L1=diag(-(F1+[B1(ngrid),B1(1:ngrid-1)])) ...
   +diag(F1(1:ngrid-1),-1)+diag(B1(1:ngrid-1),1);
L1(1,ngrid) = F1(ngrid);
L1(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(1,ngrid) = F2(ngrid);
L2(ngrid,1) = B2(ngrid);

M = zeros(2*ngrid,2*ngrid);
M(1:ngrid,1:ngrid) = L1;
M(ngrid+1:2*ngrid,ngrid+1:2*ngrid) = L2;
M = M-diag([k12,k21])...
    +diag(k21,ngrid)+diag(k12,-ngrid);

%Compute the matrix need for the Crank-Nicholson routine
G = eye(2*ngrid,2*ngrid) - dt/2 * M;
H = eye(2*ngrid,2*ngrid) + dt/2 * M;

HG = G\H;

%Integrate the probabilities
for i = 1:nsteps
    p = HG*p;
end

%compute the density
rho = p/dx;