clear all

% Analyzed data from Britten and Kohne, Science, 1968.
%
% Executes properly on MATLAB version 9.4.0.813654 (R2018a).
% Executes poorly on OCTAVE version 4.2.1.
% Joshua C. Vaughan 5/14/2020

% Get points out of a figure plot. Initially used the code, below, for
% this.

% basepath = pwd;
% img_nm_in = [basepath '\Britten_key_plot.png'];
% img = imread(img_nm_in);
% 
% figure(1); clf
% set(gcf,'position',[817 322 1092 635])
% image(img);
% 
% %[xc,yc] = ginput(4); % mark the 4 corners of plot in image
% xc_im = [74  72 886 884];
% yc_im = [476 122 121 477];
% 
% 
% %% Get some inputs. Assign variables after batches if desired.
% % Code here adapted from an online source.
% % https://www.mathworks.com/matlabcentral/answers/403814-getting-points-with-ginput-from-image-and-plotting-simultaneously-each-selected-point
% 
% n = 0;
% while true
%    [x, y, button] = ginput(1);
%    if isempty(x) || button(1) ~= 1; break; end
%    n = n+1;
%    x_n(n) = x(1); % save all points you continue getting
%    y_n(n) = y(1);
%    hold on
%    plot(x(1), y(1), 'r')
%    drawnow
% end
% 
% % copy data to clipboard using command clipboard('copy',x_n) or
% clipboard('copy',y_n) then paste within this .m file.

%%%
% points from the image for defining x and y axes, units in pixels of image
x1_im = 73;
x2_im = 885;
y1_im = 122;
y2_im = 477;

% defining points on the plot
x1_pl = -7; % log10 output
x2_pl = 4.7; % log10 output
y1_pl = 0;
y2_pl = 1;

% slope for plot units over image units
m_x = (x2_pl-x1_pl)/(x2_im-x1_im);
m_y = (y2_pl-y1_pl)/(y2_im-y1_im);

%%% Fit data with a second-order rate law 
s = fitoptions('Method','NonlinearLeastSquares',...
    'Lower',[-inf],...                                 % parameter lower bound [A_min B_min]
    'Upper',[inf],...                                 % parameter upper bound [A_max B_max]
    'StartPoint',[1]);                                % estimate of fit parameters
f = fittype('1./(1+A*x)','options',s);                 % form of the fit (exponential, here)
t = logspace(-7,5,2^8);                             % time vector for evaluating fit

%%% mouse satelite curve
x_im_msat = [297.6613  307.8049  318.8307  328.9743  335.1487  343.9693 347.4975  355.8770  359.4052  362.4924  366.0207  371.3130  379.2515  384.1028 388.9541  399.5388  407.4773  414.5337 418.0620  422.9133  427.3236];
y_im_msat = [179.5557  198.9610  242.1817  272.1716  298.6333  324.6540 334.3566  342.2951  356.8490  377.5773  382.4287  404.0390 407.5673  422.5622 430.9417  437.5572  446.3777  443.2905 451.2290  450.7880  456.9624];

% convert units
x_msat = x1_pl + m_x*(x_im_msat-x1_im);
y_msat = y1_pl + m_y*(y_im_msat-y1_im);

% fit to second order integrated rate law
x_msat_lin = 10.^x_msat;                            % convert x to linear scale
y_msat_lin = 1-y_msat;                              % convert y to denatured concentration
FitStruct_msat = fit(x_msat_lin',y_msat_lin',f);    % do the fit here
y_msat_fit = feval(FitStruct_msat,t);               % evaluate the fit

%%% MS2 bacteriophage
x_im_ms2 = [329.015957446808 357.21012783116 362.952878345916 378.37626544269 383.954937371311 398.886088709677 413.817240048044 427.271684111187 432.686277453672 447.289271619767 447.617428792038 445.812564344543 460.907794269046 469.111723575841 489.785625428963 495.200218771448 531.133429135209];
y_im_ms2 = [156.837949537422 171.62057212203 197.289044513714 205.97683516936 248.625989297081 255.339282076445 298.383335779423 332.739598826753 343.401887358683 366.700962298827 385.656141911148 402.241924071928 429.884894339895 425.540999012072 467.795253564536 453.183969280039 467.795253564536];

% convert units
x_ms2 = x1_pl + m_x*(x_im_ms2-x1_im);
y_ms2 = y1_pl + m_y*(y_im_ms2-y1_im);

% fit to second order integrated rate law
x_ms2_lin = 10.^x_ms2;                            % convert x to linear scale
y_ms2_lin = 1-y_ms2;                              % convert y to denatured concentration
FitStruct_ms2 = fit(x_ms2_lin',y_ms2_lin',f);    % do the fit here
y_ms2_fit = feval(FitStruct_ms2,t);               % evaluate the fit

%%% T4 bacteriophage
x_im_t4 = [419.947216026081 434.143295298559 453.301111015786 470.253710535347 477.145011153054 491.341090425532 498.370217055594 504.158909574468 511.877166266301 520.008900995196 524.970637439945 528.278461736445 533.653676218257 537.237152539465 541.509758922443 544.40410518188 547.436277453672 548.814537577213 554.051926046671 565.215833047358 571.004525566232 579.136260295127 590.300167295813 606.287984728895];
y_im_t4 = [129.969372041418 139.618243591063 161.142649355656 180.069282010729 198.624804221585 221.262541318829 245.755830637159 261.713579738495 276.186887062962 298.082403271772 302.906839046595 310.329047930937 324.431244811187 335.193447693484 346.697871464214 354.120080348557 360.428957900248 369.706719005676 375.273375668932 390.117793437617 401.622217208348 410.899978313775 423.14662297294 440.588813851145];

% convert units
x_t4 = x1_pl + m_x*(x_im_t4-x1_im);
y_t4 = y1_pl + m_y*(y_im_t4-y1_im);

% fit to second order integrated rate law
x_t4_lin = 10.^x_t4;                            % convert x to linear scale
y_t4_lin = 1-y_t4;                              % convert y to denatured concentration
FitStruct_t4 = fit(x_t4_lin',y_t4_lin',f);    % do the fit here
y_t4_fit = feval(FitStruct_t4,t);               % evaluate the fit

%%% e. coli
x_im_ecoli = [549.198116849691 552.5409445779 567.583669354839 577.194299073439 584.576376973233 590.287041008922 593.769153225806 602.544076012354 606.86189516129 613.686835106383 617.865369766644 624.551025223061 633.047379032258 645.58298301304 654.636474776939 661.322130233356 683.468363932738 703.386045813315 714.250235929993 724.83585706932];
y_im_ecoli = [152.174741682019 156.785075277485 174.842215193062 189.441604912039 204.040994631016 217.103606484838 226.70846814206 248.223358254237 263.591136905792 281.648276821369 294.326694208902 312.383834124479 338.893252298411 366.939448337499 387.301755050809 402.285339236075 430.331535275163 450.693841988473 455.30417558394 458.76192578054];

% convert units
x_ecoli = x1_pl + m_x*(x_im_ecoli-x1_im);
y_ecoli = y1_pl + m_y*(y_im_ecoli-y1_im);

% fit to second order integrated rate law
x_ecoli_lin = 10.^x_ecoli;                            % convert x to linear scale
y_ecoli_lin = 1-y_ecoli;                              % convert y to denatured concentration
FitStruct_ecoli = fit(x_ecoli_lin',y_ecoli_lin',f);    % do the fit here
y_ecoli_fit = feval(FitStruct_ecoli,t);               % evaluate the fit

%%% calf
x_im_calf = [661.238632463967 668.654984557309 675.082489704873 680.356340082361 685.135766986959 688.596731297186 697.166738160604 703.099819835278 703.429435483871 707.22001544269 706.725591969801 710.680979752917 714.30675188744 719.250986616335 722.052719629375 726.502530885381 732.270804735758 735.237345573095 738.863117707618 742.324082017845 744.631391557996 747.597932395333 750.234857584077 753.36620624571 755.508707961565 760.123327041867 764.737946122169 765.891600892244 770.011796499657 774.461607755662 778.911419011668 781.053920727522 800.995667467399 804.786247426218 806.928749142073 808.247211736445 814.015485586822 818.959720315717 821.431837680165 824.727994166095 826.870495881949 830.331460192176 832.638769732327 835.934926218257 840.549545298559 843.845701784489 844.504933081675 846.977050446122 848.460320864791 849.943591283459 849.119552161977 853.404555593686 854.723018188058 857.689559025395 859.832060741249 859.832060741249 860.656099862732];
y_im_calf = [122.221252061834 123.753272100782 123.370267091045 125.285292139729 122.987262081308 129.881352256572 130.264357266309 131.41337229552 125.285292139729 121.838247052097 134.094407363678 136.009432412362 132.179382314993 131.41337229552 137.158447441573 131.796377305256 136.392437422099 144.052537616837 149.031602743417 146.350567675259 152.478647831049 157.840717967365 173.543923366578 171.245893308157 174.692938395789 182.736043600264 191.928163833949 194.992203911844 205.716344184477 213.376444379215 227.54762973948 234.441719914745 286.913406248699 286.913406248699 294.956511453174 297.637546521332 311.042721862124 323.681887183441 340.917112621602 350.109232855287 359.684358098709 359.684358098709 361.982388157131 365.8124382545 371.94051841029 388.792738838713 397.601854062662 399.516879111347 401.431904160031 402.963924198979 408.325994335295 407.942989325558 427.859249831877 419.433139617665 419.050134607928 422.880184705297 431.306294919509];

% convert units
x_calf = x1_pl + m_x*(x_im_calf-x1_im);
y_calf = y1_pl + m_y*(y_im_calf-y1_im);

% fit to second order integrated rate law
x_calf_lin = 10.^x_calf;                            % convert x to linear scale
y_calf_lin = 1-y_calf;                              % convert y to denatured concentration
FitStruct_calf = fit(x_calf_lin',y_calf_lin',f);    % do the fit here
y_calf_fit = feval(FitStruct_calf,t);               % evaluate the fit

%%% polyUandA
FitStruct_poly = FitStruct_calf;
FitStruct_poly.A = 1/2e-6;
y_poly_fit = feval(FitStruct_poly,t);               % evaluate the fit

%%% Make plot that looks like the one in the paper. 
ax_lw = 1;
pl_lw = 2;
pl_ms = 20;
ax_fs = 20;
la_fs = 24;
ti_fs = 24;
fontname = 'Arial';


%%% plot a couple of decays and fits
figure(1); clf
plot(t,y_calf_fit,'c',t,y_ecoli_fit,'m',t,y_t4_fit,'b',t,y_ms2_fit,'g',t,y_msat_fit,'r','linewidth',pl_lw)
hold on
plot(x_ecoli_lin,y_ecoli_lin,'m.',x_t4_lin,y_t4_lin,'b.',x_ms2_lin,y_ms2_lin,'g.',x_msat_lin,y_msat_lin,'r.','markersize',pl_ms)
hold off
xlim([0 2])
ylim([0 1.05])
xlabel('[SS(0)]\times{\itt}  (M s)','fontsize',la_fs)
ylabel('[SS({\itt})]/[SS(0)]','fontsize',la_fs)
set(gca,'ytick',[0 1])
set(gcf,'color',[1 1 1])
set(gca,'linewidth',2)
set(gca,'fontsize',ax_fs)
set(gca,'tickdir','out')
set(gca,'xminortick','on')
% set(gca,'xtick',xt)
% set(gca,'xticklabel',xtl)
set(gca,'fontname',fontname)
h = legend('cow','{\itE. coli}','T4','MS2','mouse sat.');
set(h,'box','off')
set(h,'fontsize',20)
set(h,'position',[0.5622    0.4007    0.3033    0.3220])
set(gca,'LooseInset',get(gca,'TightInset'));
