clear all;
close all;

%time span for integration of the selection system
time=0:1200;

%defined start values in M
%start value for Target protein concentration
T0=0;
%start value for Evobody concentration
E0=0; 
%start value for DNA
D0=0.000000067; 
%start value for Evobody-Target complex
ET0=0;
%start value for DNA-Target complex
TD0=0;
%start value for Evobody-Target-DNA complex
ETD0=0;

%%
%start of the selection system

%different association rates for Evobody and Target protein
k_assET = [379000,379000,379000];
%different disassociation rates for Evobody and Target protein
k_dissET = [0.00272,0.10347,1.153]; 

%expression rate of plasmid 1 in 1/(s*M)
k_exprP1 = 7.38017 * 10^-10;
%expression rate of plasmid 2 in 1/(s*M)
k_exprP2 = 7.38017 * 10^-10;

%concentration of ampicillin in medium [M]
amp_med = 0.000572; 

%vector of all reproduction rates
growth_rates = [];


%calls  for the computations of the selection system
%loop length represents the number of different adapted E. coli cells
loop_length = length(k_assET);

allAmpPeri=[];
allbetaLac=[];

for i=1:loop_length

    figure(2*i-1);
    %Solves the differential equation
    [T,X] = ode23s(@bilbico , time, [E0 T0 D0 ET0 TD0 ETD0 0 0 0 amp_med], [],k_assET(i), k_dissET(i), k_exprP1, k_exprP2 );

    
%     plot(T,X);
%     legend('E','T','D','ET', 'TD','ETD', 'bla_Cyt', 'bla_Peri', 'amp_Peri', 'amp_Medium');
%     
%     
%     xlabel('time (sec))', 'FontSize', 15) ;
%     ylabel('M', 'FontSize', 15) ;
%     %test=legend('100% affinity', '95% affinity', '40% affinity');
%     
%     %save figure
%     print(['-f' int2str(2*i-1)],['bild' int2str(2*i-1)],'-dpng');
    figure(2*i);
    
    
    %remeber all beta-lactamse
    allbetaLac(:,i)=X(:,8);
    
    %Amp concentration in periplasm for t = 0 to t = 1200 
    amp_Peri = X(:,9);
    allAmpPeri(:,i)=amp_Peri;
    display(amp_Peri(end));
    plot(T,amp_Peri);
    growth_rates(i) = doubling_ecoli_time( amp_Peri(end) );
end
%plot beta-lactamase  concentration
% figure(300);
% i=2*i+1;
% plot(T,allbetaLac);
% 
% legend('optimal affinity','moderate affinity','bad affinity');
% xlabel('time (sec))', 'FontSize', 15) ;
% ylabel('M', 'FontSize', 15) ;
% title('beta-lactamase concentration');
% 
% %save figure
% %  print(['-f' int2str(i)],'Bla' ,'-dpng');



%plot amp peri
% figure(301);
% i=i+1;
% plot(T,allAmpPeri);
% 
% legend('optimal affinity','moderate affinity','bad affinity');
% xlabel('time (sec))', 'FontSize', 15) ;
% ylabel('M', 'FontSize', 15) ;
% 
% title('amp concentration');
% 
% save figure
%     print(['-f' int2str(i)],'Amp','-dpng');
    


%%
%growth curves
figure(300);
hold on;

cc=hsv(3);

%time span for growth curves
t = 0:20;
color=['b','g','r'];
%plots the different good adapted E. coli cultures with the calculated growing rates 
for i=1:loop_length
    display(growth_rates(i));
    %Solves the differential equation
    [T,y] = ode45(@verhulst_growth_dynamicR, t, [growth_rates(i) 0.1], [] , 0.38);
    plot(T,y(:,2),color(i));
    
    
end

legend('optimal affinity','moderate affinity','bad affinity');
xlabel('time (h))', 'FontSize', 15) ;
ylabel('OD', 'FontSize', 15) ;
title('Growth curves');

%%
%competitive growth curves

%time for  in h
t=[1,2,3,4,5,6,7,8];
% t=[1,2];

%inserte Od
insert_od=0.1;

%All start values
start_values=[];

%gives all different adapted E. coli cells the same start value
for i=1:loop_length
    start_values(i)=insert_od/loop_length;

end

%how much better the best must be
difference_percent=[10,20,30,40,50,60,70,80,90,95,100];
% difference_percent=[10,20];

%how many ⌀fferent percent values
how_many=length(difference_percent);

%all need reinoculation
all_reinoculation=[];

%save all number of reinoculations for the different parameter
reinoculation=[];

%
for r= 1:length(t)
    %test when the best is better 10%, 20%, .., 90%, 100%
    for i=1:how_many
        display(i);
        %compute how many reinoculationsare used 
        [batches,who]=find_time_when_best_is_percent_better( t(r), insert_od, start_values, growth_rates, difference_percent(i) );
        all_reinoculation(i)=batches;
        reinoculation(i,r)=batches;
    
    end

%     display(all_reinoculation);
%     display(growth_rates);

    figure(100+r);

    plot(difference_percent, all_reinoculation);

end

%plot 3d
figure(100+r+1);
mesh(   t,difference_percent, reinoculation );
xlabel('time between inoculation [h] ');
ylabel('percentage of best adapted culture');
zlabel('number of reinoculations');
%save figure
print(['-f' int2str(100+r+1)],['bild3D' int2str(r+1)],'-dpng');

%plot 3D
f=figure(100+r+2);
surf( t,difference_percent, reinoculation );
xlabel('time between inoculation [h] ','FontSize', 30);
ylabel('percentage of best adapted culture','FontSize', 30);
zlabel('number of reinoculations','FontSize', 30);
% set(get(gca,'xlabel'),'rotation',30);
% set(get(gca,'ylabel'),'rotation',-30);
%save figure
print(['-f' int2str(100+r+2)],['bild3d' int2str(r+2)],'-dpng');

