% Numerical solution of ternary FRAP model with solvent, bleached and
% unbleached species. Model is assumed to be equilibrated
% (bleached+unbleached=const.=pt). Then bleached species initial
% conditions are introduced. Integration of model via pdepe.
pa = '/Users/hubatsch/ownCloud/Dropbox_Lars_Christoph/DropletFRAP/FRAP_paper/';
a = -1;
b = 0.000025;
u0 = 0.05;
e = 0.4;
e_g0 = 0.16;
o_g0 = 0.2;
%%
g0 = gamma0(x, a, b, e_g0, o_g0);
pt = phi_tot(x, a, b, e, u0);
%% Make useful mesh (by inverting the tanh profile and using this as spacing)
x = linspace(-a-1, -a+1, 3000);
g = gamma0(x, a, 1*b, e_g0, o_g0);
g_unique = unique(g);
x = linspace(g_unique(1), g_unique(end-1), 300);
g_inv = spacing(x, a, 2*b, e_g0, o_g0);
g_inv = g_inv(2:end-1);
x = [linspace(-a-1, g_inv(2), 30), g_inv(3:end-2), ...
    linspace(g_inv(end-1), -a+1, 30), linspace(-a+1.1, 300, 3000)];
%% Solve pde
tic
t = linspace(0, 2, 200);
fh_ic = @(x) flory_ic(x, a, u0);
fh_bc = @(xl, ul, xr, ur, t) flory_bc(xl, ul, xr, ur, t, u0);
fh_pde = @(x, t, u, dudx) flory_hugg_pde(x, t, u, dudx, a, b, e, u0,...
                                         e_g0, o_g0);
sol = pdepe(0, fh_pde, fh_ic, fh_bc, x, t);
toc
%% Plotting
figure(1); hold on;
for i = 1:length(t)
    cla; xlim([-a-1, -a+3]); ylim([0, 0.7]); 
    ax = gca;
    ax.FontSize = 16;
    xlabel('position'); ylabel('volume fraction');
    plot(x, phi_tot(x, a, b, e, u0), 'LineWidth', 2, 'LineStyle', '--'); 
    plot(x, sol(i, :), 'LineWidth', 2);
    shg
    pause();
%     print([num2str(i),'.png'],'-dpng')
end
%% Boundary condition: Jump makes sense
[~, ind] = min(abs(x+a)); % Find x position of boundary
plot(max(sol(:, 1:ind-1), [], 2)./min(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(u0+e-sol(:, 1:ind-1), [], 2)./max(u0-sol(:, ind+1:end), [], 2), 'LineWidth', 2)
make_graph_pretty('time', 'ratio outside/inside', 'bleached component');
%% 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')
%% Function definitions for pde solver
function [c, f ,s] = flory_hugg_pde(x, t, u, dudx, a, b, e, u0,...
                                    e_g0, o_g0)
% Solve with full ternary model. Analytical derivatives.
% pt ... phi_tot
% gra_a ... analytical gradient of phi_tot

pt = phi_tot(x, a, b, e, u0);
gra_a = gradient_analytical(x, a, b, e);
g0 = gamma0(x, a, b, e_g0, o_g0);
c = 1;
f = g0*(1-pt)/pt*(pt*dudx-u*gra_a);
s = 0;
end

function u = flory_ic(x, a, u0)
    % FRAP initial condition
    if x < -a; u = 0.0; else; u = u0; end
    % Peak outside initial condition
%     if (x < -a+0.2) && (x > -a+0.1); u = 10; else; u = 0; end
%     u = 0.3;
%     u = phi_tot(x, -50, 1);
%%
%% Frank's solution to the transfer/rate problem via Laplace transform
x0 = 3;
D_m = 0.11;%g0(1)*(1-u0-e)/(u0+e); % to make equal to ternary FRAP
D_p = 0.342;%g0(end)*(1-u0)/u0;
ga = 1/9;
p_out = @(D_p, D_m, ga, x0, x, t) 1./(2*sqrt(D_p*pi*t))*...
            (exp(-(x+x0).^2./(4*D_p*t))*(ga*sqrt(D_p)-sqrt(D_m))./...
            (ga*sqrt(D_p)+sqrt(D_m))+exp(-(x-x0).^2./(4*D_p*t)));
if x > -a
    u = p_out(D_p, D_m, ga, x0, x, 3/100);
else
    u = 0;
end
end

function [pl,ql,pr,qr] = flory_bc(xl, ul, xr, ur, t, u0)
    pl = 0;
    ql = 1;
%     % No flux
%     pr = 0;%ur - 0.01;
%     qr = 1;
    % Dirichlet BC
    pr = ur - u0;
    qr = 0;
end

function g0 = gamma0(x, a, b, e_g0, o_g0)
    g0 = e_g0*(tanh((x+a)/b)+1)/2+o_g0;
end

function sp = spacing(x, a, b, e_g0, o_g0)
    sp = b*atanh(2/(e_g0)*(x-o_g0)-1)-a;
end

function p = phi_tot(x, a, b, e, u0)
    p = e*(tanh(-(x+a)/b)+1)/2+u0;
end

function grad = gradient_analytical(x, a, b, e)
    grad = -e*(1-tanh(-(x+a)/b).^2)/b*0.5;
end