% 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.

a = -50;
b = 0.01;
c = 1/100000;
u0 = 0.3;
%% Solve pde
x = [linspace(49, 55, 600), linspace(55.01, 200, 300)];
t = linspace(0, 0.001, 1000);
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, c, u0);
sol = pdepe(0, fh_pde, fh_ic, fh_bc, x, t);
%% Plotting
figure(1); hold on;
for i = 1:length(t)
    cla; xlim([49, 53]); ylim([0, 1.5]);
    plot(x, phi_tot(x, a, b, u0)); plot(x, sol(i, :)); pause(0.03);
end

%% Plot and check derivatives of pt
figure; hold on;
x = linspace(40, 60, 100);
plot(x, phi_tot(x, -50, 0.25));
plot(x, gradient_analytical(x, -50, 0.25));

%% Function definitions for pde solver
function [c, f ,s] = flory_hugg_pde(x, t, u, dudx, a, b, c_p, u0)  
% Solve with full ternary model. Analytical derivatives.
% pt ... phi_tot
% gra_a ... analytical gradient of phi_tot

pt = phi_tot(x, a, b, u0);
gra_a = gradient_analytical(x, a, b);
c = c_p;
f = ((u0+1.3)-pt)/pt*(pt*dudx-u*gra_a);
s = 0;
end

function u = flory_ic(x, a, u0)
    if x < -a
        u = 0.0;
    else
        u = u0;
    end
%     u0 = 0.3;
%  u0 = phi_tot(x, -50, 1);
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 p = phi_tot(x, a, b, u0)
    p = (tanh(-(x+a)/b)+1)/2+u0;
end

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