%% 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');

%% Figures
% Time course
figure; hold on;
xlim([0, 3]); ylim([0, 1]); 
ax = gca;
ax.FontSize = 16;
xlabel('x [\mum]'); ylabel('volume fraction'); 
plot(T.x, T.sol(1:1:100, :), '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, 0), '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', {-1, b(7/3, 10^-6), 0.5, e(7/3),...
                   0, 1, 100, 0, 0, 'Constituent'},...
                   linspace(0.0, 4000, 5000), 1);
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');

%%
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));
%% Check dependence of simulation accuracy on precision parameter
color = ['r', 'b', 'k', 'm', 'y', 'g'];
figure; hold on;
prec = [0.1, 0.5, 1.1, 2, 5, 10];
t = linspace(0.0, 40, 500);
tic
for i = 1:length(wi)
Ts{i} =Ternary_model(2, 'FRAP', [-1, 0.02, 0.5, 0.3, 0.3, 0.5, 40], t, prec(i));
Ts{i}.solve_tern_frap();
% 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]);


%% Can clearly see trench going away from minimum which becomes really flat
nu = 10^-5;
chi = 7/3;
b = nu^(1/3)*sqrt(chi/(chi-2));
e = sqrt(3/8*(chi-2));
b = 0.05;
e = 0.4;
e_g0 = 0.05:0.1:4.8;
u_g0 = 0.05:0.1:4.8;
T = cell(length(e_g0), length(u_g0));
tic
for jj = 1:length(u_g0)
    parfor i = 1:length(e_g0)
    T{i, jj} = Ternary_model(2, 'FRAP', [-1, b, 0.5, e, e_g0(i), u_g0(jj), 40],...
                      linspace(0.0, 100, 500), 1);
    T{i, jj}.solve_tern_frap();
    end
end
toc
%% Set reference simulation
e_ref = 3;
u_ref = 40;
dist_squ = cell(length(e_g0), length(u_g0));
figure; hold on;
for jj = 1:length(u_g0)
    for i = 1:length(e_g0)
    dist_squ{i, jj} = sum((sum(T{i, jj}.sol(:, 1:100), 2)-...
                           sum(T{e_ref, u_ref}.sol(:, 1:100), 2)).^2);
%     dist_squ{i, jj} = sum((sum(T{i, jj}.sol, 2)-sum(T{e_ref, u_ref}.sol, 2)).^2);
    plot(sum(T{i, jj}.sol, 2))
    pause()
    end
end
%%
surf(u_g0, e_g0, cell2mat(dist_squ));
xlabel('u_g0'); ylabel('e_g0');
h = gca;
set(h,'zscale','log')


%% Swapping tau_in and tau_out does not result in the same intensity curves
t = Ternary_model(2, 'FRAP', [-1, 0.001, 0.5, 0.4, 0.02, 0.55, 40],...
                      linspace(0.0, 100, 500), 1);
t.solve_tern_frap()
%% Inside and outside time scale
Xi = t.phi_t(1)/t.phi_t(end);
R_eff = Xi^(1/3) * (-t.a);
D_in = (1-t.phi_t(1))*Ternary_model.gamma0(t.x(1), t.a, t.b, t.e_g0, t.u_g0);
D_out = (1-t.phi_t(end))*Ternary_model.gamma0(t.x(end), t.a, t.b, t.e_g0, t.u_g0);
tau_out = R_eff^2/D_out;
tau_in = t.a^2/D_in;
%% Turn around tau_in and tau_out
D_in1 = tau_out/t.a^2;
gamma0_in = D_in1/(1-t.phi_t(1));
D_out1 = R_eff/tau_in;
gamma0_out = D_out1/(1-t.phi_t(end));
%%
t1 = Ternary_model(2, 'FRAP', [-1, 0.001, 0.5, 0.4, gamma0_in, gamma0_out, 40],...
                      linspace(0.0, 100, 500), 1);
t1.solve_tern_frap();
%%
figure; hold on;
plot(sum(t.sol(:, 1:100), 2))
plot(sum(t1.sol(:, 1:100), 2))
axis([0, 250, 0, inf])
%%
cost = to_min_theory(t.sol(:, 1:100), t.t', [0.02, 0.55],...
                     sum(t.sol(:, 1:100), 2));
%% Trench instead of minimum, sim doesn't show second minimum in this case
x_seed = [0.01, 0.6];
x_seed = [10, 0.01];
x_seed = [34.0290    0.1761];
f_min = @(y) to_min_theory(t.sol(:, 1:100), t.t', y,...
                     sum(t.sol(:, 1:100), 2));
opt = optimset('MaxFunEvals', 2000, 'PlotFcns',@optimplotfval);
for i = 1:1
    x_seed = fminsearch(f_min, x_seed, opt);
end
%%
e_g0 = 55:5:105;
u_g0 = 0.05:0.05:0.4;
T = cell(length(e_g0), length(u_g0));
tic
for jj = 1:length(u_g0)
    parfor i = 1:length(e_g0)
    T{i, jj} = Ternary_model(2, 'FRAP', [-1, 0.001, 0.5, 0.4, e_g0(i), u_g0(jj), 40],...
                      linspace(0.0, 100, 500), 1);
    T{i, jj}.solve_tern_frap();
    end
end
toc
%% Compare with reference simulation from above
dist_squ = cell(length(e_g0), length(u_g0));
figure; hold on;
for jj = 1:length(u_g0)
    for i = 1:length(e_g0)
    dist_squ{i, jj} = sum((sum(T{i, jj}.sol(:, 1:100), 2)-...
                           sum(t.sol(:, 1:100), 2)).^2);
%     dist_squ{i, jj} = sum((sum(T{i, jj}.sol, 2)-sum(T{e_ref, u_ref}.sol, 2)).^2);
    plot(sum(T{i, jj}.sol, 2))
%     pause()
    end
end
%%
surf(u_g0, e_g0, cell2mat(dist_squ));
xlabel('u_g0'); ylabel('e_g0');
h = gca;
set(h,'zscale','log')