%% 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];
%% Try out different precisions
prec = [0.25, 0.5, 1, 2, 3.5, 5];
parfor i = 1:length(prec)
    tic
    T_prec(i) = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, ...
                              e(7.7/3), 0, 1, 300, 7], t, prec(i));
    T_prec(i).solve_tern_frap()
    toc
end

%% Same precisions, different starting positions
x0 = [-0.05, -0.03, -0.01, 0.01, 0.03, 0.05];
parfor i = 1:length(x0)
T_diff_x0(i) = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, ...
                             e(7.7/3), 0, 1, 300, 6+x0(i)], t, 0.5);
T_diff_x0(i).solve_tern_frap()
end
%%
for i = 1:length(T)
    csvwrite(['x_X73_', num2str(x0(i)), '.csv'], T_diff_x0(i).x+T(i).a);
    csvwrite(['sol_X73_', num2str(x0(i)), '.csv'], T_diff_x0(i).sol(:, :));
end

%% Frank's/Stefano via Laplace transform vs Fokker Planck
t = linspace(0, 10, 1000);
tic
T_mov = Ternary_model(0, 'Gauss', [-6, b(7/3, 10^-15), 0.5, e(7/3),...
                   0, 1, 300, 7, 0], t, 0.2);
T_mov.solve_tern_frap()
toc
%% 
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);
%% Moving boundary
t = linspace(0, 10, 300);
tic
T_mov = Ternary_model(0, 'Uniform', [-10, b(7/3, 10^-6), 0.5, e(7/3),...
                   0, 1, 10, 7, 0.5], t, 0.02);
T_mov.solve_tern_frap();
toc
%%
csvwrite('MovingBound.csv', [T_mov.x', T_mov.sol(1, :)', T_mov.sol(end, :)'])
%% FRAP jump length
t = linspace(0, 5, 300);
T_mov = Ternary_model(0, 'FRAP', [-5, b(7/3, 10^-6), 0.5, e(7/3),...
                   0, 1, 10, 7, 0], t, 0.2);
T_mov.solve_tern_frap();
% csvwrite('FRAP.csv', [T_mov.x', T_mov.sol(1, :)', T_mov.sol(end, :)'])

%% Calculate fluxes FRAP
figure(1);
for i = 1:300
x_interp = (T_mov.x(1:end-1)+T_mov.x(2:end))/2;
u = T_mov.sol(i, :);
u_interp = (u(1:end-1)+u(2:end))/2;
pt = Ternary_model.phi_tot(x_interp, T_mov.a, T_mov.b, T_mov.e,...
                           T_mov.u0, T_mov.v*T_mov.t(i));
gra_a = Ternary_model.gradient_analytical(x_interp, T_mov.a, T_mov.b,...
                                          T_mov.e, T_mov.v*T_mov.t(i));
g0 = Ternary_model.gamma0(x_interp, T_mov.a+T_mov.v*T_mov.t(i), T_mov.b,...
                          T_mov.e_g0, T_mov.u_g0, T_mov.v*T_mov.t(i));
dudx = diff(u)./diff(T_mov.x);
hold on; cla;
plot(x_interp, g0.*(1-pt)./pt.*(pt.*dudx-u_interp.*gra_a))
pause()
end
% csvwrite('FRAP.csv', [T_mov.x', T_mov.sol(1, :)', T_mov.sol(end, :)'])
%% Plot with full ternary model
T_temp = T_mov;
D_m = 1*(1-T_temp.u0-T_temp.e); % to make equal to ternary FRAP
D_p = 1*(1-T_temp.u0+T_temp.e);
ga = (T_temp.u0-T_temp.e)/ (T_temp.u0+T_temp.e);
for i = 1:100%length(T_mov.t)
    cla;
    figure(1); hold on;
    xlim([-2, 2]);
    ylim([0, 1]);
    plot(T_temp.x+T_temp.a, T_temp.sol(i, :), 'r.');%, 'LineWidth', 2);
    plot(x_left, p_in(D_p, D_m, ga, T_temp.x0+T_temp.a, x_left,...
        t(i)+0.01/1000), 'k-.', 'LineWidth', 2);
    plot(x_right, p_out(D_p, D_m, ga, T_mov.x0+T_temp.a, x_right,...
        t(i)+0.01/1000), 'k-.', 'LineWidth', 2)
    make_graph_pretty(['x [' char(956) 'm]'], 'c [a.u.]', '')
%     print([num2str(i),'.png'],'-dpng')
    shg; pause();
end
%% Check integral of solution, should be mass conserving and sum to 1.
% integrate to 50, to avoid right boundary
[~, ind] = min(abs(T_prec(4).x-50));
bins = (T_prec(4).sol(:, 1:ind-1)+T_prec(4).sol(:, 2:ind))/2;
bin_size = diff(T_prec(4).x(1:ind));
sum(bins(4, :).*bin_size)
%%
figure; hold on;
T_prec(1).plot_sim('plot', 10, 'magenta')
T_prec(2).plot_sim('plot', 10, 'red');
T_prec(3).plot_sim('plot', 10, 'green');
T_prec(4).plot_sim('plot', 10, 'k');
T_prec(5).plot_sim('plot', 10, 'blue');
%% Check whether partitioning factor is kept throughout simulation.
% This should work for steep boundaries.
for i = 1:4
pks_min = findpeaks(-T_prec(5).sol(i, :));
pks_max = findpeaks(T_prec(5).sol(i, :));
pks_max(1)/pks_min(1)
end

%% Solve integrals
a = 5;
nu = 10^-6;
chi = 7.7/3;
b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
e = @(chi) sqrt(3/8*(chi-2));
t = [0, 0.05, 0.1];
% T1 = Ternary_model(0, 'Gauss', [-5, b(chi, nu), 0.5, e(chi), ...
%                      0, 1, 10, a], t, 0.2);
x0 = 5.0:0.002:7.01;
%%
T = {};
parfor i = 1:length(x0)
tic
b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
e = @(chi) sqrt(3/8*(chi-2));
T{i} = Ternary_model(0, 'Gauss', [-a, b(chi, nu), 0.5, e(chi), ...
                     0, 1, 10, x0(i)], t, 0.2);
T{i}.solve_tern_frap();
toc
end
% save prob_laplace_X_7_7_short_2
%%
ls = 0.0005:0.005:2;
p = nan(1, length(ls));
parfor i = 1:length(ls)
p(i) = int_prob(ls(i), T, x0, 0.1);
end
%% Do we need to look at the left side as well?
figure; hold on;
plot(ls, p/N);
% plot(lx, xt);
% plot(ls, q/50000);
%%
for i = 1:2:length(T)
    T{i}.plot_sim('plot', 3, 'red')
end
%%
int_prob(0.3, T, x0, 0.1)
%%
int_prob_simple(0.01, T, x0)
%% Normalization factor
tic
N = normalization(T, x0,0.1);
toc
%% Write to file
csvwrite('jump_length_7_7_lb01.csv', ls)
csvwrite('prob_7_7_lb01.csv', p/N);
%% distribution for x0 can be taken from phi_tot (steady state)
function p = int_prob(l, T, x0, lb)
delta_x0 = diff(x0);
p = 0;
for i = 1:length(delta_x0)
    x = (x0(i)+x0(i+1))/2;
    if x-5>lb
        if x-l > 5-lb; break; end
        p_i = @(j) interp1(T{j}.x, T{j}.sol(3, :), x-l);
        p2 = (p_i(i)+p_i(i+1))/2;
        p = p + delta_x0(i)*...
                T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0)*p2;
    end
end
end

function p = int_prob_simple(l, T, x0)
p = 0;
p_i = @(j, x0) interp1(T{j}.x, T{j}.sol(3, :), x0-l);
for i = 1:length(x0)
    if x0(i)-l > 5; break; end % Is this right?
    % Calculate left bin boundary
    if i == 1
        left_bound = -T{1}.a; 
    else
        left_bound = (x0(i)-x0(i-1))/2;
    end
    % Calculate right bin boundary
    if i == length(x0)
        right_bound = T{1}.system_size;
    else
        right_bound = (x0(i+1)+x0(i))/2;
    end
    bin_width = right_bound-left_bound;
    
    p = p + bin_width*...
            T{1}.phi_tot(x0(i), T{1}.a, T{1}.b, T{1}.e, T{1}.u0)*...
            p_i(i, x0(i));
end
end

function p = normalization(T, x0, lb)
delta_x0 = diff(x0);
logi = T{1}.x<5-lb;
delta_x = diff(T{1}.x(logi));
p_i = zeros(1, length(delta_x0));
for i = 1:length(delta_x0)
    sol = T{i}.sol(3, logi);
    sol = (sol(1:end-1)+sol(2:end))/2;
    x = (x0(i)+x0(i+1))/2;
    if x > 5+lb
        p_x0i = T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0);
        p_i(i) = delta_x0(i)*sum(sol.*delta_x)*p_x0i;
    end
end
p = sum(p_i);
end