% 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, 40, 100), linspace(40.4, 60, 50), linspace(60.4, 100, 100)];
x = linspace(40, 60, 300);
% x = linspace(0.5, 100.5, 1001);

t = linspace(0, 500000, 50);
% y = linspace(0.3000001, 1.2999999, 500);
% x = 0.5*atanh(2*y-1.6)+50;
a = -50;
b = 0.5;
% x_t = linspace(0.5, 101.5, 1011);
% global phi_tot_g
% phi_tot_g = phi_tot(x_t, b, a);
% global der_phi 
% der_phi = gradient_analytical(x_t, a, b);
% global lapl_phi
% lapl_phi = laplacian_analytical(x_t, a, b);

% 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_a, @flory_ic, @flory_bc, x, t);

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

%% Plot and check derivatives of pt
figure; hold on;
% x = linspace(40, 60, 100);
plot(x, phi_tot(x, -50, 0.5));% plot(x, gra_pt(x, -50, 0.5, 0.001)); plot(x, lap_pt(x, -50, 0.5, 0.001));
% plot(x, gralap_pt(x, -50, 0.5, 0.001)); plot(x, laplap_pt(x, -50, 0.5, 0.001));
plot(x, gradient_analytical(x, -50, 0.5));
plot(x, laplacian_analytical(x, -50, 0.5));
%%
phi_tot(x, -50, 1)
%% Function definitions for pde solver
function [c, f ,s] = flory_hugg_a(x, t, u, dudx)  
% Solve with full ternary model. Analytical derivatives.
% pt ... phi_tot
% gra_a ... analytical gradient of phi_tot
% lap_a ... analytical laplacian of phi_tot

pt = @(x) phi_tot(x, -50, 1);
gra_a = @(x) gradient_analytical(x, -50, 1);
lap_a = @(x) laplacian_analytical(x, -50, 1);
% gra_a = @(x) gra_pt(x, -50, 0.5, 0.0001);
% lap_a = @(x) lap_pt(x, -50, 0.5, 0.0001);
c = 1/(1.3-pt(x));
f = dudx;
s = u/(1.3-pt(x))*(lap_a(x)+(gra_a(x)/pt(x))^2-lap_a(x)/pt(x))...
    -dudx/(1.3-pt(x))/pt(x)*gra_a(x);
end

% function [c, f ,s] = flory_hugg_a(x, t, u, dudx)  
% % Solve with full ternary model. 
% global phi_tot_g
% global der_phi 
% global lapl_phi
% 
% i = round(10*x);
% % disp(i)
% c = 1/(1-phi_tot_g(i));
% f = dudx;
% s = u/(1-phi_tot_g(i))*(lapl_phi(i)+(der_phi(i)/phi_tot_g(i))^2-lapl_phi(i)/phi_tot_g(i))...
%     -dudx/(1-phi_tot_g(i))/phi_tot_g(i)*der_phi(i);
% end

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

function gpt = gra_pt(x, a, b, delta)
    gpt = (phi_tot(x+delta, a, b)-...
           phi_tot(x-delta, a, b))/(2*delta);
end

function lpt = lap_pt(x, a, b, delta)
    lpt = (gra_pt(x+delta, a, b, delta)-...
           gra_pt(x-delta, a, b, delta))/(2*delta);
end
% 
% function glpt = gralap_pt(x, a, b, delta)
%     glpt = (lap_pt(x+delta, a, b, delta)-...
%             lap_pt(x-delta, a, b, delta))/(2*delta);
% end
% 
% function llpt = laplap_pt(x, a, b, delta)
%     llpt = (gralap_pt(x+delta, a, b, delta)-...
%             gralap_pt(x-delta, a, b, delta))/(2*delta);
% end

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

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