% Calculate RBS strengths from GFP expression data
function rbs_calc2() 

% matlab code exported by Gen5Parser
Experiment1_26_7_10;

strains = {
    'K398500A'
    'K398501A'
    'K398502A'
    'K398503A'
    'K398504A'
    'K173000A'
    'I13401'
    'LB+Amp'
};

rbsname = {
    'J61100'
    'J61101'
    'J61107'
    'J61117'
    'J61127'
    'B0032'
    'fluor blank'
    'abs blank'
};

% plot all graphs of J61100
egfprange = 1:96;
odrange = 97:192;
egfpod = data(:, egfprange) ./ data(:, odrange);

% Normalize by subtracting LB+Amp OD curves from OD data
% Plot LB+Amp
% plot (x, data(:, H1_OD600_600_index:H12_OD600_600_index));
for i=97:12:192
    data(:, i:i+11) = data(:, i:i+11) - data(:, H1_OD600_600_index:H12_OD600_600_index);
end
% Same for EGFP
for i=1:12:96
    data(:, i:i+11) = data(:, i:i+11) - data(:, H1_EGFP_485_20_528_20_index:H12_EGFP_485_20_528_20_index);
end


figure(1);
x = (1:size(data,1)).* ( 1/6 );
plot (x, egfpod(:, 1 : 12*8));
title('EGFP/OD');

figure(2);
plot (x, data(:, 96 + (1 : 12*8)));
title('OD');

figure(3);
plot (x, data(:, 1:12 + 12*7));
title('EGFP');

figure(5);
plot (x, data(:, H1_EGFP_485_20_528_20_index:H12_EGFP_485_20_528_20_index));
title('blank EGFP');

% plot mean and std.dev. for each RBS
figure(4);
%[output_mean, output_std] = plotexp(data, 1:12:72);
%legend(rbsname);

% .5 hour till 3 hour: 
measure_range = 5:15;
alpha_deg = log(2)/(41*60);
for i = 1:72
    timestep = (measure_range(end)-measure_range(1)) * 10;
    od1 = data(measure_range(1), 96+i);
    od2 = data(measure_range(end), 96+i);
    
    gfp1 = data(measure_range(1), i);
    gfp2 = data(measure_range(end), i);
    
    mu = 1 / timestep * (log(od2) - log(od1));
    dt = log(2) / mu; % doubling time

    alpha = 1/mu + alpha_deg;
    beta = .5 * alpha * (gfp2^2 - gfp1^2) / (gfp2 - gfp1);
    
    fprintf('GFP production rate/min: %f\n', beta);
    prod(i) = beta;
    %fprintf('Growth: %f; %f\n', mu, dt);
end

mprod_std = mean(prod((1:12) + (6-1)*12));
for i=1:6
    mprod = mean(prod((1:12) + (i-1)*12));
    
    strength = .3 * (mprod / mprod_std);
    fprintf('%s RBS strength = %f\n', rbsname{i}, strength);
end


end
