%% Solve Fokker Planck equation for given parameter set
% Define parameters for tanh, according to Weber, Zwicker, Julicher, Lee.
b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
e = @(chi) sqrt(3/8*(chi-2));
t = [0, 0.1, 1, 9.9];
%%
T = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
                  0.2, 300, 7], t, 0.5);
T.solve_tern_frap()
%%
T1 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
                   0.2, 300, 7], t, 1);
T1.solve_tern_frap()
%%
T2 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
                   0.2, 300, 7], t, 2);
T2.solve_tern_frap()
%%
tic
T3 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
                   0.2, 300, 7], t, 5);
T3.solve_tern_frap()
toc
%%
T4 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
                   0.2, 300, 7], t, 10);
T4.solve_tern_frap()
%% different x0
x0 = [-0.05, -0.03, -0.01, 0.01, 0.03, 0.05];
for i = 1:length(x0)
T(i) = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
                   0.2, 300, 6+x0(i)], t, 1);
T(i).solve_tern_frap()
end
%%
for i = 1:length(T)
    csvwrite(['x_X73_', num2str(x0(i)), '.csv'], T(i).x+T(i).a);
    csvwrite(['sol_X73_', num2str(x0(i)), '.csv'], T(i).sol(:, :));
end
%% Frank's solution to the transfer/rate problem via Laplace transform
j = 4;
x0 = .01;
D_m = 1*(1-T(j).u0-T(j).e); % to make equal to ternary FRAP
D_p = 1*(1-T(j).u0+T(j).e);
ga = (T(j).u0-T(j).e)/ (T(j).u0+T(j).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)))*2;
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))*2;
x_left = linspace(-4, 0, 1000);
x_right = linspace(0, 4, 1000);
%% Plot with full ternary model
for i = 2:length(T(j).t)  
    figure(1); hold on; %cla;
    plot(x_left, p_in(D_p, D_m, ga, x0, x_left, t(i)));
    plot(x_right, p_out(D_p, D_m, ga, x0, x_right, t(i)));
    pause()
    plot(T(j).x+T(j).a, 2*T(j).sol(i, :), 'LineWidth', 2);
    axis([-2, 2, 0, 5.8]);
    shg; pause();
end
%% Check integral of solution, should be mass conserving and sum to 1.
[~, ind] = min(abs(T3.x-50)); % integrate up to 50, to avoid right boundary
bins = (T3.sol(:, 1:ind-1)+T3.sol(:, 2:ind))/2;
bin_size = diff(T3.x(1:ind));
sum(bins(4, :).*bin_size)
%%
figure; hold on;
T.plot_sim('plot', 1, 'blue');
T1.plot_sim('plot', 1, 'magenta')
T2.plot_sim('plot', 1, 'red');
T3.plot_sim('plot', 1, 'green');
T4.plot_sim('plot', 1, 'k');
%% Check whether partitioning factor is kept throughout simulation.
% This should work for steep boundaries.
for i = 1:4
pks_min = findpeaks(-T4.sol(i, :));
pks_max = findpeaks(T4.sol(i, :));
pks_max(1)/pks_min(1)
end