clear; close all;
percentc60 = 50;%Percentage viable hopping sites in the mixed layer
ebarrier = 100;%Energy barrier to entry of the mixed layer from the neat fullerene layer in meV
kt = 26;%RT thermal energy in meV
d=100;%Lattice diemsions in number of lattice sites (assumed 1 nm per site)
l=200;
w=7;
randomblend = rand(d,l,w); %gerate a random matrix with dimension representing depth length width, respectively
randomblend(randomblend<percentc60/100)=1; %Convert the correct percentage to 1 (indicates C60 sites)
randomblend(randomblend<1)=0;%rest set to zero (indicates non-fullerene sites)
randomblend(1:10,:,:)=1;%set all the neat layer to fullerene (top ten lattice sites)

N=10000;%Number of initial charges
charges = ones(N,3);%Initial charge position array (all start at lattice site (1,1,1), which is at the top of the neat layer
collected = zeros(N,1);%preallocate arrays to keep track of which charges are collected or quenched 
quenched = zeros(N,1);
totalcollected = 0;
totalquenched = 0;
skip = 10;%Number of steps to take before plotting
step = 1;
numinneatlayer = N;
while totalcollected < N/2 %Keep Looping until half of the charges are collected (for speed, the percentage collected when finishing doesn't affect statistical result)
    for j = 1:N %loop over all charges
        if collected(j)==0 %if not collected
            stepdir = 6*rand;%choose a random number between 0,6
            if stepdir < 1 %0-1 indicates hopping attempt down
                if charges(j,1)==d-1 %If at quenching interface, quench
                    quenched(j)=1;
                    totalquenched = totalquenched + 1;
                    collected(j) = 1;
                    totalcollected = totalcollected + 1;
                elseif randomblend(charges(j,1)+1,charges(j,2),charges(j,3))~=0 %if the site below is a viable hopping site (not non-fullerene)
                    if charges(j,1) == 10 && rand()<exp(-ebarrier/kt) %if at the neat/mixed interface, hopping succeeds with boltzmann probability
                        charges(j,1) = charges(j,1)+1;
                    elseif charges(j,1) ~= 10 %otherwise hop down with 100% chance of success
                        charges(j,1) = charges(j,1)+1;
                    end
                end
            elseif stepdir < 2 %2-6 represent the remaining directions. No barrier to hopping except at boundaries
                if charges(j,1)>1 && randomblend(charges(j,1)-1,charges(j,2),charges(j,3))~=0
                    charges(j,1) = charges(j,1)-1;
                end
            elseif stepdir < 3
                if charges(j,2)==l-1
                    collected(j)=1;
                    totalcollected = totalcollected + 1;
                elseif randomblend(charges(j,1),charges(j,2)+1,charges(j,3))~=0
                    charges(j,2) = charges(j,2)+1;
                end
            elseif stepdir < 4
                if charges(j,2)>1 && randomblend(charges(j,1),charges(j,2)-1,charges(j,3))~=0
                    charges(j,2) = charges(j,2)-1;
                end
            elseif stepdir < 5
                if charges(j,3) < w && randomblend(charges(j,1),charges(j,2),charges(j,3)+1)~=0
                    charges(j,3) = charges(j,3)+1;
                elseif charges(j,3) == w && randomblend(charges(j,1),charges(j,2),1)~=0
                    charges(j,3)=1;
                end
             elseif stepdir < 6
                if charges(j,3) > 1 && randomblend(charges(j,1),charges(j,2),charges(j,3)-1)~=0
                    charges(j,3) = charges(j,3)-1;
                elseif charges(j,3) == 1 && randomblend(charges(j,1),charges(j,2),w)~=0
                    charges(j,3)=w;
                end
            end
        end
    end
    if mod(step,skip)==0 %periodically plot to visualize charge distribution
        figure(1)
        plot3(charges(:,1),charges(:,2),charges(:,3),'*')
        axis([1 d 1 l 1 w])
        title(['Collection Efficiency = ' num2str((totalcollected-totalquenched)/totalcollected) ', P_{up} = ' num2str(ebarrier) ', ' num2str(percentc60) '% C60, % Finished = ' num2str(totalcollected/(N/2))])
        view(2)
        drawnow
    end
    step = step + 1;
    numinneatlayer(step)=sum(charges(:,1)<11);
end
efficiency = (totalcollected-sum(quenched))/totalcollected  %calculate efficiency of charge collection