%% 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;
plot(s(1:10:400, :)', 'Color', [0.12156   0.4666666   0.7058823], 'LineWidth', 1.5);
plot(v_fit(1:10:400, :)', '--', 'Color', [1.   0.49803   0.0549], 'LineWidth', 1.5);
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, s(:, end))
make_graph_pretty('t [s]', 'intensity [a.u.]', '');
print([T.pa, 'Timecourse_BC'],'-depsc')
%% 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', 'ratio outside/inside', 'unbleached component');
axis([0, 200, 0, 0.035])