%concentration can only be plotted for two dimensions, in this case x and
%y, as the third is reserved for concentration itself
%D = 4.9×10^-6 cm2 s^-1 from

%http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=1797063
%assuming bacteria population doubles every 30 minutes
%returns the required population size in 1 microlitre in the form of a
%sphere. 

%this piece, as you can see takes the inputs from the user
clear all
time2 = input('time over which you want to look at HSL concentration (minutes) =');
scale = input('size of bacteria population increment per iteration =');
start = input('number at which to begin iteration =');
step1 = input('number of bacteria population iterations =');

%this sets out all the parameters
steps = time2*60;          %converting time to seconds
volume = 0.001^3;          %volume of 1 millimetre cubed in metres cubed (microlitre)
r = ((3*volume)/(4*pi))^(1/3);
area = 4*pi*(r^2);        % surface area of a cube
D = 4.9*(10^-10);          %m^2s^-1 at unknown temp
out = (D*area)/(volume*r);                   %this calculates what I assume to be the coefficient for HSL that can diffuse per second out the cube
N = zeros(steps,(step1-start));
flux = zeros(steps,(step1-start));
conc = zeros(steps,(step1-start));
start1=start+1;

gashvar = (455*volume)/(5.65*10^-17);
j = 0;

for n = start1:step1
for t = 2:steps
    N(t,(n-start)) = 49.91*n*scale + N(t-1,(n-start)) - flux(t-1,(n-start));
    flux(t,(n-start)) = N(t,(n-start))*out;
    conc(t,(n-start)) = N(t,(n-start))/((6.02214179*10^23)*volume);
    if N(t,(n-start))>=gashvar
        j = j + 1;
        quorum_time(j) = t;
        population(j) = n*scale;
        break
    end 
end

end





graph = input('Do you want a graph? (0 = no, 1 = yes)');
if graph == 1
    
n = (start:(step1-1))*scale;
t = 1:steps;


figure
xlabel('bacteria number')
ylabel('time (seconds)')
zlabel('concentration')

HSL = 455/((5.65*10^-17)*(6.02214179*10^23))
Percent_Volume = (min(population)*(5.65*10^-17)*100)/volume

mesh(n,t,conc)
hold on
plot3(min(population),max(quorum_time),HSL,'+g')
xlabel('bacteria number')

figure
plot(population,quorum_time);
ylabel('response time')
xlabel('bacteria number')
end


    



%49.1 molecules of HSL are produced per cell

%I apologise if this is confusing, feel free to e-mail me on: rory.macdonald.06@aberdeen.ac.uk if you have any questions.