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

%% Solve pde
x = linspace(0,500, 6000);
t = linspace(0, 500, 50);
% opt = odeset('RelTol',1e-14, 'AbsTol', 1e-16,'MaxStep',1e-2); 
% sol = pdepe(0,@flory_pde, @flory_ic, @flory_bc,x,t, opt);
sol = pdepe(0,@flory_hugg_pde, @flory_ic, @flory_bc,x,t);

%% Plotting
figure(1); hold on;
for i = 1:50
    cla; plot(x, pt(x)); plot(x, sol(i, :)); pause(0.04);
end

%% Plot and check derivatives of pt
figure; hold on;
x = linspace(140, 170, 500);
plot(x, pt(x)); plot(x, gra_pt(x, 0.001)); plot(x, lap_pt(x, 0.001));
plot(x, gralap_pt(x, 0.001)); plot(x, laplap_pt(x, 0.001));
plot(x, gra_a(x, -155, 0.5));
plot(x, lap_a(x, -155, 0.5));
% pt(0)
% gra_pt(0, 0.00001)
% lap_pt(0, 0.00001)
% lap_a(0, -155, 0.5)
% gralap_pt(0, 0.00001)
% laplap_pt(0, 0.00001)
%% Function definitions for pde solver
function [c, f ,s] = flory_hugg_pde(x, t, u, dudx)  
% Solve with full ternary model.
X = 3.2;
K = 1.5;
c = 1/(pt(x)+1);
f = dudx;
s = u*(lap_pt(x, 0.0001)-2*X*(lap_pt(x, 0.0001)-...
        gra_pt(x, 0.0001)^2-pt(x)*lap_pt(x, 0.0001))+...
        K*(-laplap_pt(x, 0.0001)+...
            gra_pt(x, 0.0001)*gralap_pt(x, 0.0001)+...
            pt(x)*laplap_pt(x, 0.0001)))/(pt(x)+1)+...
    dudx*(-2*X*(gra_pt(x, 0.0001)-pt(x)*gra_pt(x, 0.0001))-...
          K*gralap_pt(x, 0.0001)+K*pt(x)*gralap_pt(x, 0.0001))/...
          (pt(x)+1);
end

function u0 = flory_ic(x)
    if x<5
        u0 = 0.0;
    else
        u0 = 0.5;
    end
%     u0 = pt(x)*0.5;
end

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

function [c, f ,s] = flory_pde(x, t, u, dudx) 
% Tentatively solve system where binary mixture is assumed, with variable
% '1' (now called pt). Basically pt is fixed and now phi_u and
% phi_b can swap places but nothing else.
c = 1/pt(x);
f = dudx;
s = -u*lap_pt(x)/pt(x);
end

function p = pt(x)
    p = (tanh(-(x-155)*2)+1)/2+0.3;
end

function gpt = gra_pt(x, delta)
    gpt = (pt(x+delta)-pt(x-delta))/(2*delta);
end

function lpt = lap_pt(x, delta)
    lpt = (gra_pt(x+delta, delta)-...
           gra_pt(x-delta, delta))/(2*delta);
end

function glpt = gralap_pt(x, delta)
    glpt = (lap_pt(x+delta, delta)-...
            lap_pt(x-delta, delta))/(2*delta);
end

function llpt = laplap_pt(x, delta)
    llpt = (gralap_pt(x+delta, delta)-...
            gralap_pt(x-delta, delta))/(2*delta);
end

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

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