ternary_frap.m

%% Boundary condition: Jump makes sense
[~, ind] = min(abs(T.x+T.a)); % Find x position of boundary
plot(max(T.sol(:, 1:ind-1), [], 2)./min(T.sol(:, ind+1:end), [], 2), 'LineWidth', 2)
make_graph_pretty('time', 'ratio outside/inside', 'unbleached component');
figure; % Shouldn't be symmetric, different SS for phi_u and phi_b??
plot(min(T.u0+T.e-T.sol(:, 1:ind-1), [], 2)./max(T.u0-T.sol(:, ind+1:end), [], 2),...
    'LineWidth', 2)
make_graph_pretty('time', 'ratio outside/inside', 'bleached component');
%% Plot flux at boundary

%% Figures
% Time course
figure; hold on;
xlim([0, 3]); ylim([0, 0.5]); 
ax = gca;
ax.FontSize = 16;
xlabel('x [\mum]'); ylabel('volume fraction'); 
plot(T.x, T.sol(1:2:200, :), 'LineWidth', 2, 'Color', [135/255  204/255 250/255]);
plot(T.x, Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0), 'LineWidth', 4,...
    'LineStyle', '--', 'Color', [247/255, 139/255, 7/255]);
print([T.pa, 'ternary_time_course'], '-depsc');
%% Plot and check derivatives of pt
figure; hold on;
x = linspace(40, 60, 100000);
plot(T.x, Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0));
plot(T.x, Ternary_model.gradient_analytical(T.x, T.a, T.b, T.e));
plot(T.x(1:end-1)+mean(diff(T.x))/2, ...
     diff(Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0)/mean(diff(T.x))));
plot(T.x, Ternary_model.gamma0(T.x, T.a, T.b, T.e_g0, T.o_g0));
axis([0, 3, 0, 1])
figure;
plot(Ternary_model.gamma0(T.x, T.a, T.b, T.e_g0, T.o_g0),...
     Ternary_model.spacing(Ternary_model.gamma0(T.x, T.a, T.b, T.e_g0, T.o_g0),...
     T.a, T.b, T.e_g0, T.o_g0));
%% Fit with experimental BCs
T = Ternary_model(2, 'FRAP', [], linspace(0.0, 400, 5000));
T.solve_tern_frap();
x_seed = 60;
s = T.sol(1:end, 1:30)./max(T.sol(1:end, 1:30))/1.022;
ts = T.t(1:end)';
f_min = @(x) to_minimize(s, ts, x, 'simple_drop', s(:, end));
opt = optimset('MaxFunEvals', 2000, 'PlotFcns',@optimplotfval);
x_seed = fminsearch(f_min, x_seed, opt);
%%
[cost, v_fit, r_n, v_fit_end] = to_minimize(s, ts, x_seed, 'simple_drop', s(:, end));
%% Plot for 
figure; hold on;
s_all = T.sol(1:end, 1:400)./max(max(T.sol(1:end, 1:30)))/1.022;
plot(T.x(1:340), s_all(1:10:200, 1:340), 'Color', [0.5294, 0.8, 0.98], 'LineWidth', 2.5);
plot(T.x(1:30), v_fit(1:10:200, :)', '.', 'Color', [0.83, 0.23, 0.09], 'MarkerSize', 10);
make_graph_pretty('$x [\mu m]$', 'intensity [a.u.]', '');
% print([T.pa, 'Timecourse_model_only'],'-depsc')
%% Plot boundary condition (read out)
figure; hold on;
plot(ts(1:200), s(1:200, end), 'Color', [0.1608    0.5412    0.2118])
make_graph_pretty('t [s]', 'intensity [a.u.]', '');
print([T.pa, 'Timecourse_BC'],'-depsc')
%% Fit boundary condition with exponential
f = fit(ts(1:200), s(1:200, end), 'B+C-B*exp(-x/tau)-C*exp(-x/tau2)');
plot(f)
%% Plot flux at boundary (x close to 1) Only works close to boundary!
% T = Ternary_model(2, 'FRAP', [], linspace(0.0, 400, 2000));
% T.e_g0 = 0.001;
% T.b = 0.001;
% T.solve_tern_frap();
[~, ind_in] = min(abs(T.x+T.a+0.02)); % Find x position of boundary close to 0.9
% Calculate local diffusion coefficient
g0 = Ternary_model.gamma0((T.x(ind_in)+T.x(ind_in+1))/2, T.a, T.b, T.e_g0, T.o_g0);
pt = Ternary_model.phi_tot((T.x(ind_in)+T.x(ind_in+1))/2, T.a, T.b, T.e, T.u0);
D_in = g0*(1-pt);
Diffs_in = D_in*diff(T.sol')/(T.x(ind_in+1)-T.x(ind_in));
[~, ind_out] = min(abs(T.x+T.a-0.02)); % Find x position of boundary close to 0.9
% Calculate local diffusion coefficient
g0 = Ternary_model.gamma0((T.x(ind_out)+T.x(ind_out+1))/2, T.a, T.b, T.e_g0, T.o_g0);
pt = Ternary_model.phi_tot((T.x(ind_out)+T.x(ind_out+1))/2, T.a, T.b, T.e, T.u0);
D_out = g0*(1-pt);
Diffs_out = D_out*diff(T.sol')/(T.x(ind_out+1)-T.x(ind_out));
figure; hold on;
plot(Diffs_in(ind_in, :))
plot(Diffs_out(ind_out, :));
make_graph_pretty('time', '\bf{j}', 'boundary flux inside/outside of drop');
axis([0, 200, 0, 0.035])
%% Check out inner and outer length scales
e_g0 = [0.01, 0.05, 0.1, 0.2, 0.5, 3, 10, 100];
for i = 1:8
%     T1{i} = Ternary_model(2, 'FRAP', [], linspace(0.0, 400, 5000));
%     T1{i}.e_g0 = e_g0(i);
%     T1{i}.solve_tern_frap();
    T1{i}.calc_time_scales();
end
%%
figure; hold on;
for i = 1:length(T1)
    [~, ind] = min(abs(T1{i}.x+T1{i}.a)); 
    plot(T1{i}.t, sum(T1{i}.sol(:, 1:ind), 2)); axis([0, 100, 0, 80]);
end
%% Mix
o_g0 = [0.01, 0.05, 0.1, 0.2, 0.5, 3, 10, 100];
for i = 1:8
    T2{i} = Ternary_model(2, 'FRAP', [], linspace(0.0, 400, 5000));
    T2{i}.o_g0 = o_g0(i);
    T2{i}.solve_tern_frap();
    T2{i}.calc_time_scales();
end
%%
figure; hold on;
for i = 1:length(T1)
    [~, ind] = min(abs(T2{i}.x+T2{i}.a)); 
    plot(T2{i}.t, sum(T2{i}.sol(:, 1:ind), 2)); axis([0, 100, 0, 80]);
end
%%
%%
o_g0 = [0.01, 0.05, 0.1, 0.2, 0.5, 3, 10, 100];
for i = 1:8
    T3{i} = Ternary_model(2, 'FRAP', [], linspace(0.0, 400, 5000));
    T3{i}.e_g0 = 0.01;
    T3{i}.o_g0 = o_g0(i);
    T3{i}.solve_tern_frap();
    T3{i}.calc_time_scales();
end
%%
figure; hold on;
for i = 1:length(T1)
    [~, ind] = min(abs(T2{i}.x+T2{i}.a)); 
    plot(T2{i}.t, sum(T2{i}.sol(:, 1:ind), 2)); axis([0, 100, 0, 80]);
end
%% Outside time scale
fol = '~/ownCloud/Dropbox_Lars_Christoph/DropletFRAP/Latex/Figures/SpatialAnalysisAdvantages/';
p = logspace(1, 3, 20);
R_in = 1;
tau_out = p.^(1/3);
plot(p, tau_out);
make_graph_pretty('partitioning coefficient \xi', '\tau_{out}', '');
print([fol, 'partitioning_vs_tau'],'-depsc');
%%
T_large_p = Ternary_model(2, 'FRAP', [-1, 0.025, 0.0005, 0.5-0.0005, 5, 0.16],...
                            linspace(0.0, 400, 5000));
T_large_p.solve_tern_frap()
%%
T_middle_p = Ternary_model(2, 'FRAP', [-1, 0.025, 0.0015, 0.5-0.0015, 5, 0.16],...
                            linspace(0.0, 400, 5000));
T_middle_p.solve_tern_frap()
%%
T_small_p = Ternary_model(2, 'FRAP', [-1, 0.025, 0.04, 0.5-0.04, 5, 0.16],...
                            linspace(0.0, 400, 5000));
T_small_p.solve_tern_frap()
%%
T_middle_p.plot_sim('plot', 10);
print([fol, 'middle_partitioning'],'-depsc');
T_large_p.plot_sim('plot', 10)
print([fol, 'large_partitioning'],'-depsc');
T_small_p.plot_sim('plot', 10)
print([fol, 'small_partitioning'],'-depsc');
%% Partitioning time scales outside/inside
Din_Dout=logspace(-3, 1, 100);
Xi = logspace(0, 4, 100);
[X,Y]=meshgrid(Xi,Din_Dout);
Z=Y.*X.^(2/3);
figure; hold on;
hAx=subplot(1, 1, 1);
surf(X,Y,log(Z));
set(hAx,'xscale','log');
set(hAx,'yscale','log');
colorbar
colormap(c)
shading interp
view(2)
plot(Xi, Xi.^(-2/3), 'k');
ylim([10^-3, 10]);
make_graph_pretty('Partitioning \xi', 'D_{in}/D_{out}', '');
print([fol, 'partitioning_DinDout'],'-depsc');