%% 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-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([49, 53]); ylim([0, 0.5]); 
ax = gca;
ax.FontSize = 16;
xlabel('position'); ylabel('volume fraction'); 
plot(x, sol(1:2:300, :), 'LineWidth', 2, 'Color', [135/255  204/255 250/255]);
plot(x, phi_tot(x, a, b, e, u0), 'LineWidth', 4, 'LineStyle', '--', 'Color',...
    [247/255, 139/255, 7/255]);
% print([pa, 'ternary_time_course'], '-depsc');
%% Plot and check derivatives of pt
figure; hold on;
x = linspace(40, 60, 100000);
plot(x, phi_tot(x, a, b, e, u0));
plot(x, gradient_analytical(x, a, b, e));
plot(x(1:end-1)+mean(diff(x))/2, ...
     diff(phi_tot(x, a, b, e, u0)/mean(diff(x))));
plot(x, gamma0(x, a, b, e_g0, o_g0));
figure;
plot(gamma0(x, a, b, e_g0, o_g0),...
     spacing(gamma0(x, a, b, e_g0, o_g0), a, b, e_g0, o_g0));
%% Fit with experimental BCs
s = sol(1:end, 1:30)./max(sol(1:end, 1:30))/1.022;
ts = 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:2:100, :)', 'Color', [0.12156   0.4666666   0.7058823], 'LineWidth', 1.5);
% plot(v_fit(1:2:100, :)', '--', 'Color', [1.   0.49803   0.0549], 'LineWidth', 1.5);
make_graph_pretty('$x [\mu m]$', 'intensity [a.u.]', '');
print('Timecourse_model_only.png','-dpng')