%% 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
start = [141, 202, 240];
stop = [242, 139, 18];
c = [linspace(start(1), stop(1), 256)', linspace(start(2), stop(2), 256)', linspace(start(3), stop(3), 256)']/256;
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');
h = colorbar('northoutside', 'XTickLabel',{'10^{-6}','10^{-4}','10^{-2}','1','10^2','10^4','10^6','10^8'}, ...
            'YTick',-6:2:8);
ylabel(h, [char(964) '_{in}/' char(964) '_{out}']);
% set(h,'YTick',[-6:2:8], '10);
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');
%%
nu = 10^-15;
chi = 7/3;
b = nu^(1/3)*sqrt(chi/(chi-2));
e = sqrt(3/8*(chi-2));
t = linspace(0.0, 40, 500);
T = Ternary_model(2, 'FRAP', [-1, b, 0.5, e, 1, 0.5], t);
tic
T.solve_tern_frap();
toc
%%
x_q = linspace(0, 0.9, 100);
v_q = interp1(T.x, T.sol', x_q);
v_q = v_q/max(v_q(:));
% plot(sum(v_q, 1))
x_fit = [10, 2];
model = 'simple_drop';
f_min = @(x_fit) intensity_fit(v_q', t', x_fit, model, sum(v_q, 1));
opt = optimset('MaxFunEvals', 2000, 'PlotFcns',@optimplotfval);
for i = 1:1
    x_fit = fminsearch(f_min, x_fit, opt);
end
%%
x_fit = [8, 2];
intensity_fit(v_q', t', x_fit, model, sum(v_q, 1))
[~, v_fit] = to_minimize(v_q', t', x_fit, model);
figure; hold on;
plot(v_q);
plot(v_fit')
figure; hold on;
plot(sum(v_q));
plot(sum(v_fit, 2));
%%
color = ['r', 'b', 'k', 'm', 'y', 'g'];
figure; hold on;
wi = [0.1, 0.5, 1.1, 5, 10, 2];
tic
for i = 1:length(wi)
Ts{i} =Ternary_model(2, 'FRAP', [-1, 0.02, 0.5, 0.3, 0.3, 0.5], t);
Ts{i}.create_mesh(wi(i));
Ts{i}.solve_tern_frap();
% T.create_mesh(i);
% plot(T.x, color(i), 'Marker', '.');
% length(T.x)
toc
end
toc
axis([0, 1000, 0, 2])
%%
figure; hold on;
Ts{1}.plot_sim('plot', 10, 'red')
Ts{2}.plot_sim('plot', 10, 'blue')
Ts{3}.plot_sim('plot', 10, 'black')
Ts{4}.plot_sim('plot', 10, 'yellow')
Ts{5}.plot_sim('plot', 10, 'magenta')
Ts{end}.plot_sim('plot', 10, 'green')
xlim([0, 2]); ylim([0, 1]);