%% Solve Fokker Planck equation for given parameter set
% nu = 10^-6;
% chi = 8.1/3;
% Cubic equation solution: https://en.wikipedia.org/wiki/Cubic_equation
clear all; clc;
e = 0.7;
b = 0.025;
a1 = 1;
b1 = -8/3*e^2*b^6;
c1 = 32/3;
d1 = -32/3*e^2;
del0 = b1^2-3*a1*c1;
del1 = 2*b1^3-9*a1*b1*c1+27*a1^2*d1;
C = nthroot((del1+sqrt(del1^2-4*del0^3))/2, 3);
Xi = (-1+sqrt(-3))/2;
Chi = @(k) 1/(3*a1)*(b1+Xi^k*C+del0/(Xi^k*C));
Chi(1)
Chi(2)
Chi(3)
b^3*(Chi(3)/(Chi(3)-2))^(-3/2)
%%
clear all; clc;
nu = 10^-6;
chi = 7/3;
b = nu^(1/3)*sqrt(chi/(chi-2));
e = sqrt(3/8*(chi-2));
T = Ternary_model(0, 'Gauss', [-6, b, 0.5, e, 0.16, 0.2],...
                  linspace(0.0, 400, 5000));
% T.t = 0:0.01:20;
T.t = [0, 0.1, 1, 9.9];
T.solve_tern_frap()
%% Frank's solution to the transfer/rate problem via Laplace transform
x0 = 2;
D_m = 1*(1-T.u0-T.e); % to make equal to ternary FRAP
D_p = 1*(1-T.u0+T.e);
ga = (T.u0-T.e)/ (T.u0+T.e);
%%
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)));
p_in = @(D_p, D_m, ga, x0, x, t) 1./(sqrt(pi*t)*(sqrt(D_m)+ga*sqrt(D_p)))*...
            exp(-(x-x0*sqrt(D_m/D_p)).^2/(4*D_m*t));
x_left = linspace(-4, 0, 1000);
x_right = linspace(0, 4, 1000);
%% Plot with full ternary model
for i = 1:200   
    figure(2); hold on; cla;
    j = i+2;
    plot(x_left, p_in(D_p, D_m, ga, x0, x_left, j/100));
    plot(x_right, p_out(D_p, D_m, ga, x0, x_right, j/100));
    plot(T.x+T.a, T.sol(i, :), 'LineWidth', 2);
    axis([-6, 3, 0, 2.5]);
    shg; pause();
end
%%
figure(3); hold on; cla;

plot(x_left, p_in(D_p, D_m, ga, x0, x_left, j/100), 'LineWidth', 2);
plot(x_right, p_out(D_p, D_m, ga, x0, x_right, j/100), 'LineWidth', 2);
plot(T.x+T.a, T.sol(i, :), 'LineWidth', 2);
axis([-1, 3, 0, 0.7]);