% This file is licensed under the Allen Institute Terms of Use, which can be found at:
%http://www.alleninstitute.org/Media/policies/terms_of_use_content.html

function fp_mat = generate_fp_mat(v,exp_val);

%This function generates the transition matrix to implement the 
%drift part of the Fokker-Planck equation. 
%Inputs:
%v = extended uniform discretization vector
%exp_val = exponential shift derived from characteristic equation from
%the transport equation portion: dp/dt = d(Lvp)/dv
%
%Outputs:
%fp_mat = transition matrix that implements the drift portion of the
%Fokker-Planck equation


vr = v(1);
nv = length(v);
fp_mat = eye(nv,nv);
vleft=v;
vright=[v(2:end) v(end)+1];
dvs=vright-vleft;
new_vs = v/exp_val;         %shifted voltage values after leak (drift)
qprob=1;
f_prob=1;

%re-distribute probabilities 
for binno=1:length(new_vs)-1
    leftedge=new_vs(binno);
    rightedge=new_vs(binno+1);
    newbinwidth=rightedge-leftedge;
    if (leftedge>=vr)                                          %no reversal potential crossing
        withinbins=find(vleft<=leftedge & vright>=rightedge);  %shifted bin falls entirely within an existing bin
        fullbins=find(vleft>=leftedge & vright<=rightedge);    %shifted bin entirely encompasses an existing bin
        fullbins=setdiff(fullbins,withinbins);
        leftoverhang=find(vleft<leftedge & vright>leftedge & vright<rightedge);
        rightoverhang=find(vleft<rightedge & vright>rightedge & vleft>leftedge);
        if ~isempty(withinbins) 
            fp_mat(binno,withinbins)=fp_mat(binno,withinbins)+qprob*f_prob; 
        end;
        if ~isempty(fullbins) 
            fp_mat(binno,fullbins)=fp_mat(binno,fullbins)+qprob*f_prob*dvs(fullbins)/newbinwidth; 
        end;
        if ~isempty(leftoverhang) 
            fp_mat(binno,leftoverhang)=fp_mat(binno,leftoverhang)+qprob*f_prob*(vright(leftoverhang)-leftedge)/newbinwidth; 
        end;
        if ~isempty(rightoverhang) 
            fp_mat(binno,rightoverhang)=fp_mat(binno,rightoverhang)+qprob*f_prob*(rightedge-vleft(rightoverhang))/newbinwidth; 
        end;
    else                                                      %crossing the left (negative) boundary
        if (rightedge<=vr)                                    %entire bin crosses left boundary
            fp_mat(binno,1)=fp_mat(binno,1)+qprob*f_prob;
        else
            fullbins=find(vright<=rightedge);
            rightoverhang=find(vleft<rightedge & vright>rightedge & vleft>leftedge);
            rightoverhang=setdiff(rightoverhang,fullbins);
            if ~isempty(fullbins) 
                fp_mat(binno,fullbins)=fp_mat(binno,fullbins)+qprob*f_prob*dvs(fullbins)/newbinwidth; 
            end;
            if ~isempty(rightoverhang) 
                fp_mat(binno,rightoverhang)=fp_mat(binno,rightoverhang)+qprob*f_prob*(rightedge-vleft(rightoverhang))/newbinwidth; 
            end;
            fp_mat(binno,1)=fp_mat(binno,1)+qprob*f_prob*(max(vleft(1)-leftedge,0))/newbinwidth;
        end;
    end
    fp_mat(binno,binno)=fp_mat(binno,binno)-qprob*f_prob;
end;
