%% 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.0001, 0.02, 0.05, 0.1, 1, 50];
t = linspace(0.0, 0.001, 11);
% t=0.01;
% t = logspace(-4, -2, 20);
% Folders
next = '/home/hubatsch/Nextcloud/';
mac_next = '/Users/hubatsch/Nextcloud/';
%% Comparison with Python
T1 = Ternary_model(0, 'FRAP', {-0.1, b(7/3, 10^-12), 0.5, 0.4,...
                                     0, 1, 300, 7, 0, 'Constituent', 0},...
                                     t, 1);
T1.solve_tern_frap();
csvwrite('ML_9_1.csv', [T1.x; T1.sol]);


T2 = Ternary_model(0, 'FRAP', {-0.1, b(7/3, 10^-12), 0.9, 0,...
                               -0.895, 1, 300, 7, 0, 'Constituent', 0},...
                               t, 1);
T2.x = T1.x;
T2.solve_tern_frap();
csvwrite('ML_1_9.csv', [T2.x; T2.sol]);
%%

T1.plot_sim('plot', 1, 'red');
T2.plot_sim('plot', 1, 'blue');
shg
%% Comparison with Stefano's special case Gamma=Dm/Dp analytical solution
t = linspace(0, 10  , 3);
s = csvread('~/Desktop/langevin_droplet/python_3D_Stefano.csv');
%%
T_s = Ternary_model(2, 'Gauss', {-2, b(7/3, 10^-18), 0.5, 0.353553,...
                                     0, 1, 300, 3, 0, 'Constituent', 0},...
                                     t, 5);
T_s.solve_tern_frap();
%%
T_s.plot_sim('plot', 1, 'red');
plot(s(:, 1), 9*s(:, 2)./s(:, 1).^2, 'b', 'LineWidth', 3);
axis([0, 4, 0, 1.5])
%%
T_t = T1;
ti=4;
g01 = Ternary_model.gamma0(T_t.x(657), T_t.a, T_t.b, T_t.e_g0, T_t.u_g0,...
                           0, T_t.e, T_t.u0);
pt1 = Ternary_model.phi_tot(T_t.x(656), T_t.a, T_t.b, T_t.e, T_t.u0, 0);
g02 = Ternary_model.gamma0(T_t.x(657), T_t.a, T_t.b, T_t.e_g0, T_t.u_g0,...
                           0, T_t.e, T_t.u0);
pt2 = Ternary_model.phi_tot(T_t.x(656), T_t.a, T_t.b, T_t.e, T_t.u0, 0);

g0 = (g01+g02)/2;
pt = (pt1+pt2)/2;
g0.*(1-pt)*(T_t.sol(ti, 657)-T_t.sol(ti, 656))/(T_t.x(657)-T_t.x(656))
%% Try out different precisions
% prec = [0.5, 0.5, 1, 2, 3.5, 5];
fac = [1, 5, 10];%, 10, 100];
for i = 1:length(fac)
    b = 2e-4;
    u0 = 0.5;
    P=1;
    D_i = 2;
    D_o = 20;
    a = -1;
    [e, e_g, u_g] = calc_tanh_params(fac(i)*P, D_i, fac(i)*D_o, a, b, u0);
    T_prec(i) = Ternary_model(2, 'FRAP', {-1, b, u0, ...
                              e, e_g, u_g, 300, 7, 0, 'Constituent'},...
                              t, 2);
%     T_prec(i) = Ternary_model(2, 'FRAP', {-1, b(7.7/3, 10^-5), u0, ...
%                               e(7.7/3), e_g, u_g, 300, 7, 0,...
%                               'Constituent'}, t, 0.2);
    T_prec(i).solve_tern_frap()
end

%% Test partitioning == 1 compared to natural no mobility case.
% Tg1: g0 = 1 everywhere. 
prec = linspace(0.5, 3, 6);
for i = 1:1%length(prec)
Tpt1(i) = Ternary_model(2, 'FRAP', {-1, b(7/3, 10^-12), 0.5, e(7/3),...
                                0, 1, 300, 7, 0, 'Constituent'},...
                     t, prec(i));
Tpt1(i).solve_tern_frap();
end
%%
% Rescale mobility, such that phi_tot=1 everywhere achieves same flux
Gi = 1;%(1-Tpt1.phi_t(1))/(1-Tpt1.phi_t(1));
% Go = Tpt1.phi_t(1)/Tpt1.phi_t(end)*(1-Tpt1.phi_t(end));
Go = Tpt1.phi_t(end)/Tpt1.phi_t(1)*2;
%%
for i = 1:1%length(prec)
    Tga1(i) = Ternary_model(2, 'FRAP', {-1, b(7/3, 10^-12), ...
                                        Tpt1(i).phi_t(1), 0, -0.83, Gi,...
                                        300, 7, 0,'Constituent'},...
                            t, prec(i));
    Tga1(i).x = Tpt1(i).x;
    Tga1(i).solve_tern_frap();
end
%%
for i = 1:1%length(prec)
Tpt1(i).plot_sim('plot', 1, 'g', 1, i)
Tga1(i).plot_sim('plot', 1, 'k', 1.0, i)
end
%% Calculate flux for Tpt1 at t=0
i = 6; ti=3;
fac = 1;%0.585;
[dudx, x, sol, flux, g0, pt] = calc_flux(Tpt1(i), ti);
[dudx1, x1, sol1, flux1, g01, pt1] = calc_flux(Tga1(i), ti);
figure; hold on; title('flux');
plot(x, flux); 
plot(x1, fac*flux1);
axis([0.90, 1.11, 0, inf]);
% figure; hold on; title('flux stretched');
% plot(fac*(x-1), flux); 
% plot(x1-1, fac*flux1);
% axis([-0.1, 0.11, 0, inf]);
% % figure; hold on; title('dudx');
% % plot(x, dudx); 
% % plot(x1, dudx1);
% % axis([0.99, 1.01, 0, inf]);
figure; hold on; title('Diffusion coefficients');
plot(x, g0.*(1-pt));
plot(x1, g01.*(1-pt1));
axis([0.99, 1.01, 0, inf]);
%%
figure; hold on; title('sol');
plot(x, sol); 
plot(x1, sol1);
axis([0.90, 1.11, 0, inf]);
%%
figure(3); hold on; 
% fact = [0.9965, 0.94, 1.01, 1.98];fact(i)*
c = {'g', 'm', 'k', 'r'};
for i = 1:3%length(T_prec)
%     T_prec(i).plot_sim('plot', 1, c{i})
% Is this the right normalization? %/T_prec(i).phi_t(1)
    plot(T_prec(i).x, T_prec(i).sol'/T_prec(i).phi_t(1), c{i});
    axis([0, 20, -inf, 1]);
end
axis([0, 2, 0, 1])
%% 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),...
                   19, 1, 300, 7, 0, 0, 'Constituent', 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);
%% 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 = 20*(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);
    plot(x_left, p_in(10*D_p, D_m, ga/1, T_temp.x0+T_temp.a, x_left,...
        t(i)+0.01/1000), 'b-.', 'LineWidth', 2);
    make_graph_pretty(['x [' char(956) 'm]'], 'c [a.u.]', '',...
                      [T_temp.a, 10, 0, inf])
%     print([num2str(i),'.png'],'-dpng')
    shg; pause();
end

%% %%%%%%%%%%%%%%%%% STEADY STATE JUMP LENGTH DISTRIBUTION %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
params = {-5, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0, 1, 10, 7, 0, 0,...
          'Constituent', 0};
t = [0, 0.05, 0.1, 1];
direc = 1;
x0 = sort(5-direc*(0:0.001:4.01));
%% Run simulations for 'delta' IC across outside
T = {};
parfor i = 1:length(x0)
    tic
    T{i} = Ternary_model(0, 'Gauss', params, t, 0.2);
    T{i}.x0 = x0(i);
    T{i}.solve_tern_frap();
    toc
end
% save prob_laplace_X_7_7_short_2
%% Calculate probabilities for each jump length in ls.
l_max = 4;
t_ind = 3;
ls = -direc*(0.000:0.001:l_max);
p = nan(1, length(ls));
tic
parfor i = 1:length(ls)
    tic
    p(i) = int_prob(ls(i), T, x0, direc, t_ind, 5, 0, 0.02);
    toc
end
toc
%% Normalization factor
% N = normalization(T, x0, 0, t_ind, direc, -params{1})
N=nansum(p)/length(p)*l_max;
% N = sum(p)*0.001;
m = sum(ls.*p/N)/length(p)*l_max;
% should sum to one
sum(p/N)/length(p)*l_max
%%
sum((T{1}.sol(1, 1:end-1)+T{1}.sol(1, 2:end))/2.*diff(T{1}.x))
%% Do we need to look at the left side as well?
% figure; hold on;
plot(ls, p/N);
%% Write to file
% cd([next_mac, 'Langevin_vs_MeanField/Data_Figs_FokkPla/jump_length/']);
csvwrite('jump_length_7_7_lb002.csv', ls)
csvwrite('prob_7_7_lb002.csv', p);

%% Now for d = 0.5 across the entire domain for different boundary position
params = {-5, b(7.7/3, 10^-6), 0.5, e(7.7/3), 'const', 1, 10, 7, 0, 0,...
          'Constituent', 0};
t = [0, 0.05, 0.1];
direc = 1;
bp = [3, 5, 7];

for j = 1:length(bp)
    x0 = sort(bp(j)-direc*(0:0.001:2.02));
    % Run simulations for 'delta' IC across outside
    T = {};
    parfor i = 1:length(x0)
        tic
        T{i} = Ternary_model(0, 'Gauss', [params(1:7), {x0(i)},...
                                          params(9:end)], t, 0.2);
        T{i}.solve_tern_frap();
        toc
    end
    % Calculate probabilities for each jump length in ls.
    l_max = 2;
    t_ind = 3;
    ls = -direc*(0.000:0.005:l_max);
    p = nan(1, length(ls));
    tic
    parfor i = 1:length(ls)
        tic
        p(i) = int_prob(ls(i), T, x0, direc, t_ind, bp(j), 0, 0);
        toc
    end
    toc
    save(['prob_7_7_bp', num2str(bp(j))]);
    pa = '/Users/hubatsch/Nextcloud/Langevin_vs_MeanField/data/constD/';
    pa = '';
    csvwrite([pa, 'jump_length_7_7_bp', num2str(bp(j)), '.csv'], ls)
    csvwrite([pa, 'prob_7_7_bp', num2str(bp(j)), '.csv'], p);
end

%% Now for d = 0.5 across the entire domain bp=0, lb=0.02
params = {-5, b(7.7/3, 10^-6), 0.5, e(7.7/3), 'const', 1, 10, 7, 0, 0,...
          'Constituent', 0};
t = [0, 0.05, 0.1];
direc = 1;
bp = [3, 5, 7];

for j = 1:length(bp)
    x0 = sort(bp(j)-direc*(0:0.001:3.02));
    % Run simulations for 'delta' IC across outside
    T = {};
    parfor i = 1:length(x0)
        tic
        T{i} = Ternary_model(0, 'Gauss', [params(1:7), {x0(i)},...
                                          params(9:end)], t, 0.2);
        T{i}.solve_tern_frap();
        toc
    end
    % Calculate probabilities for each jump length in ls.
    l_max = 3;
    t_ind = 3;
    ls = -direc*(0.000:0.005:l_max);
    p = nan(1, length(ls));
    tic
    parfor i = 1:length(ls)
        tic
        p(i) = int_prob(ls(i), T, x0, direc, t_ind, bp(j), 0, 0.02);
        toc
    end
    toc
    save(['prob_7_7_lb002_bp', num2str(bp(j))]);
    pa = '/Users/hubatsch/Nextcloud/Langevin_vs_MeanField/data/constD/';
    pa = '';
    csvwrite([pa, 'jump_length_7_7_lb002_bp_march', num2str(bp(j)), '.csv'], ls)
    csvwrite([pa, 'prob_7_7_lb002_bp_march', num2str(bp(j)), '.csv'], p);
end
%%
p_in = 0.853553390593274;
p_out = 1-p_in;
ratios = [1, 1.75, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30]; %
a=(1-ratios)./(p_out-ratios*p_in);
for i = 1:length(a)
    if i > 5
    [D_ratio_01(i), m_01(i)] = run_jump_lengths(-5, a(i), 7/3, 1, 4, ...
                         ['73_01_', num2str(ratios(i))], 10^-6, 3, 0.5, 1);
    end
end

for i = 1:length(a)
    if i > 5
    [D_ratio_10(i), m_10(i)] = run_jump_lengths(-5, a(i), 7/3, 1, 4, ...
                         ['73_10_', num2str(ratios(i))], 10^-6, 4, 0.5, 1);
    end
end

p_in = 0.960977222864644;
p_out = 1-p_in;
a = (1-ratios)./(p_out-ratios*p_in);

for i = 1:length(a)
    if i > 5
    [D_ratio_01_77(i), m_01_77(i)] = run_jump_lengths(-5, a(i), 7.7/3, ...
                            1, 4, ['773_01_', num2str(ratios(i))],...
                            10^-6, 3, 0.5, 1);
    end
end

for i = 1:length(a)
    if i > 5
    [D_ratio_10_77(i), m_10_77(i)] = run_jump_lengths(-5, a(i), 7.7/3, ...
                            1, 4, ['773_01_', num2str(ratios(i))],...
                            10^-6, 4, 0.5, 1);
    end
end
%% Constant flux time course
t = linspace(0, 5, 300);
T_mov = Ternary_model(0, 'Gauss', {-5, b(7/3, 10^-6), 0.5, e(7/3),...
                   0, 1, 10, 7, 0, 1, 'Const_flux', 0}, t, 0.2);
T_mov.solve_tern_frap();
%% %%%%%%%%%%%%%%%%%%% MOVING BOUNDARY TIME COURSE %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% t = linspace(0, 10, 9000); % Entropy
t = linspace(0, 1, 1001); % Jump size distro
params_mov_bound = {-10, b(7/3, 10^-6), 0.5, e(7/3),...
                    0, 1, 20, 7, 0.5, 'Client', 0};
T_mov = Ternary_model(0, 'phi_tot', params_mov_bound, t, 0);
T_mov.solve_tern_frap();
norm_fac = 1/sum(diff(T_mov.x).*...
                 (T_mov.sol(1, 1:end-1)+T_mov.sol(1, 2:end))/2);
s = '~/Nextcloud/Langevin_vs_MeanField/Data_Figs_FokkPla/';
% csvwrite([s, 'MovingBound.csv'], [T_mov.x', norm_fac*T_mov.sol(1, :)',...
%                                   norm_fac*T_mov.sol(end, :)'])
%% Flux and entropy change for moving boundary
chi_phi = -4.530864768482371;
x_interp = (T_mov.x(1:end-1)+T_mov.x(2:end))/2;
g0 = 0.5;
s_dot = zeros(1, length(T_mov.t));
for i = 1:length(T_mov.t)
    u = norm_fac*T_mov.sol(i, :);
    u_interp = (u(1:end-1)+u(2:end))/2;
    gra_a = Ternary_model.gradient_analytical(x_interp, T_mov.a, T_mov.b,...
                                              T_mov.e, T_mov.v*T_mov.t(i));
    dudx = diff(u)./diff(T_mov.x);
%     g0*(dudx+chi_phi*u*gra_a)
    f = g0*(dudx+chi_phi.*u_interp.*gra_a);
    s_dot(i) = sum(diff(T_mov.x).*f.^2./(g0.*u_interp));
    % figure(1); hold on;
    % plot(0:10, zeros(1, 11), 'LineWidth', 2);
    % plot(x_interp, f, 'LineWidth', 2)
    % set(gca,'fontsize', 18)
    % make_graph_pretty(['x [' char(956) 'm]'], 'flux', '',...
    %                   [0, T_mov.system_size, min(f(:)), max(f)])
    % pause()
end
% csvwrite([s, 'Mov_Bou_Flux.csv'], [x_interp; f])
% csvwrite([s, 'Mov_Bou_Entr.csv'], [T_mov.t; s_dot])
%%
col = 'b';
figure(2); hold on;
% plot(u, col)
% plot(u_interp, col)
% plot(gra_a, col)
plot(dudx, col)
% plot(f, col)
%% MOVING BOUNDARY JUMP LENGTH DISTRIBUTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load entropy_mov_bound_in_to_out.mat
direc = -1; % 1: jump from left to right, -1: jump from right to left.
x0 = sort(10-direc*(0:0.005:4.01));
% t_snap = repmat((0.25:0.5:9.75)', 1, length(x0)); % midpoint
t_snap = repmat((0.:0.5:9.5)', 1, length(x0)); % starting point
si_t = size(t_snap);
x0 = sort(10-T_mov.v*t_snap-direc*repmat(0:0.005:4.01, [si_t(1), 1]), 2);
%% Run simulations for 'delta' IC across outside
F = {};
t = linspace(0, 1, 21);
for j = 1:si_t(1)
    tic
    parfor i = 1:si_t(2)
        params = params_mov_bound;
        params{8} = x0(j, i);
        params{1} = -10+params{9}*t_snap(j);
        F{i, j} = Ternary_model(0, 'Gauss', params, t, 0.1);
        F{i, j}.solve_tern_frap();
    end
    toc
end

% save entropy_mov_bound_out_to_in.mat
% save entropy_mov_bound_in_to_out.mat
%%  Calc. probs. for each jump length in ls and sum over time
ls = -direc*(-0.3:0.01:2.5);
n_T=20;
p = nan(length(ls), n_T);
bound = -T_mov.a-T_mov.v*T_mov.t;
for j = 1:n_T
    tic
    F1 = F(:, j);
    bp = bound(j);
    parfor i = 1:length(ls)
        p(i, j) = int_prob(ls(i), F1, x0(j, :), direc, j, bp, 10, T_mov);
    end
    toc
end
next = '/home/hubatsch/Nextcloud/';
% csvwrite([next, 'jump_length_mov_bound_in_out.csv'], [-ls', sum(p, 2)])
% csvwrite([next, 'jump_length_mov_bound_out_in.csv', [ls', sum(p, 2)])
%% %%%%%%%%%%%%%%%%%%%%%% FRAP TIME COURSE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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, 'Constituent', 0}, t, 0.2);
T_mov.solve_tern_frap();
norm_fac = 1/sum(diff(T_mov.x)...
                 .*(T_mov.sol(1, 1:end-1)+T_mov.sol(1, 2:end))/2);
s = '~/Nextcloud/Langevin_vs_MeanField/Data_Figs_FokkPla/';
% csvwrite([s, 'FRAP.csv'], [T_mov.x', norm_fac*T_mov.sol(1, :)',...
%          norm_fac*T_mov.sol(end, :)'])
%% Flux and entropy change for FRAP
s_dot = zeros(1, length(T_mov.t));
for i = 1:length(T_mov.t)
x_interp = (T_mov.x(1:end-1)+T_mov.x(2:end))/2;
u = norm_fac*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),...
                          T_mov.e, T_mov.u0);
dudx = diff(u)./diff(T_mov.x);
f = -g0.*(1-pt)./pt.*(pt.*dudx-u_interp.*gra_a);
s_dot(i) = sum(diff(T_mov.x).*f.^2./(g0.*(1-pt).*u_interp));
% figure(1); hold on; cla;
% plot(x_interp, f)
% pause()
end
% csvwrite([s, 'FRAP_Flux.csv'], [x_interp; f])
% csvwrite([s, 'FRAP_entropy.csv'], [T_mov.t; s_dot]);
%% 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

%% %%%%%%%%%%%%%%%%%%%%%% FRAP JUMP LENGTH DISTRIBUTION %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
params = {-5, b(7/3, 10^-6), 0.5, e(7/3), 0, 1, 10, 7, 0, 'Constituent', 0};
t = linspace(0, 2, 101);
bp = 5;
direc = -2; % 1: jump from left to right, -1: jump from right to left.

%% 
[ls_i_i, p_i_i] = frap_distro(t, params, 2, bp, 1);
[ls_o_o, p_o_o] = frap_distro(t, params, -2, bp, 1);
[ls_i_o, p_i_o] = frap_distro(t, params, 1, bp, 1);
[ls_o_i, p_o_i] = frap_distro(t, params, -1, bp, 1);
%%
p_i_i(isnan(p_i_i)) = 0;
p_o_o(isnan(p_o_o)) = 0;
p_i_o(isnan(p_i_o)) = 0;
p_o_i(isnan(p_o_i)) = 0;
%%
i = 2;
figure; hold on;
plot(ls_i_i, p_i_i(:, i));
plot(ls_o_o, p_o_o(:, i));
plot(ls, p_i_o(:, i));
plot(ls, p_o_i(:, i));
%%
figure; hold on;
plot(ls_i_i, sum(p_i_i, 2));
plot(ls_o_o, sum(p_o_o, 2));
plot(ls_i_o, sum(p_i_o, 2));
plot(ls_o_i, sum(p_o_i, 2));
% save('/data/biophys/hubatsch/MatlabWorkspaces/prob_laplace_190121_2.mat');
%% Save time points for comparison with python
csvwrite('lb_in_out.csv', ls)
csvwrite('prob_in_out.csv', p);
%%
integrand = @(p, q) abs(p.*log(p./q));
integrate = @(x, y) nansum(diff(x)'.*(y(2:end)+y(1:end-1))/2);
KL = zeros(1, 50);
for i = 1:42
i_t = 1+i;
j_t = 2+i;
n_i = @(ii, oo, io, oi) abs(integrate(ls_i_i, ii) + integrate(ls_o_o, oo)+...
                            integrate(fliplr(ls_i_o), flipud(io)) +...
                            integrate(ls_o_i, oi));
N_i = n_i(p_i_i(:, i_t), p_o_o(:, i_t), p_i_o(:, i_t), p_o_i(:, i_t));
N_j = n_i(p_i_i(:, j_t), p_o_o(:, j_t), p_i_o(:, j_t), p_o_i(:, j_t));
i_i = integrate(ls_i_i, integrand(p_i_i(:, i_t)/N_i,...
                flipud(p_i_i(:, j_t))/N_j));
o_o = integrate(ls_o_o, integrand(p_o_o(:, i_t)/N_i,...
                flipud(p_o_o(:, j_t))/N_j));
i_o = integrate(fliplr(ls_i_o), integrand(flipud(p_i_o(:, i_t))/...
                N_i, p_o_i(:, j_t)/N_j));
o_i = integrate(ls_o_i, integrand(p_o_i(:, i_t)/N_i, p_i_o(:, j_t)/N_j));
KL(i) = i_i + o_o + i_o + o_i;
end

plot(t(2:51), KL/0.02);
% ylim([0, 1])
%%
figure; hold on;
plot(ls, p_i_i(:, j_t)/N_j, 'b');
plot(ls, p_o_o(:, j_t)/N_j, 'b');
plot(ls, p_i_o(:, j_t)/N_j, 'b');
plot(ls, p_o_i(:, j_t)/N_j, 'b');
%%
plot(T_mov.t, s_dot)
hold on; xlim([0, 2]);
plot(t(2:101), 0.617*abs(KL)/0.02);
%%
plot(t(2:51), 0.557*abs(KL)/0.04);
%% %%%%%%%%%%%%%%%%%%% INTEGRATION FUNCTION DEFINITIONS %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ls, p] = frap_distro(t, params, direc, bp, ind_delta)
x0 = sort(5-abs(direc)/direc*(0:0.002:4.01));

% Run simulations for 'delta' IC across outside
F = {};
parfor i = 1:length(x0)
    tic
    F{i} = Ternary_model(0, 'Gauss', params, t, 0.1);
    F{i}.x0 = x0(i);
    F{i}.solve_tern_frap();
    toc
end

%  Calc. probs. for each jump length in ls and sum over time
T_mov = Ternary_model(0, 'FRAP', params, t, 0.2);
T_mov.solve_tern_frap();
ls = -abs(direc)/direc*(0.001:0.04:4);
if abs(direc) == 2 % jumps within same phase
    ls = sort([ls, -ls]); % need pos. and neg. l, since direction not set
end
n_T=45;
p = nan(length(ls), n_T);
for j = 1:n_T
    tic
    parfor i = 1:length(ls)
        p(i, j) = int_prob(ls(i), F, x0, direc, j, bp, ind_delta, T_mov);
    end
    toc
end
% save prob_laplace_X_7_3_FRAP_in_out
end


function p = normalization(T, x0, lb, ind_t, direc, bp)
delta_x0 = diff(x0);
logi = -direc*(T{1}.x) < -direc*(bp-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(ind_t, logi);
    sol = (sol(1:end-1)+sol(2:end))/2;
    x = (x0(i)+x0(i+1))/2;
    if -direc*x > -direc*(bp+lb)
        p_x0i = T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0, 0);
        p_i(i) = delta_x0(i)*sum(sol.*delta_x)*p_x0i;
    end
end
p = sum(p_i);
end


function [dudx, x, sol, flux, g0, pt] = calc_flux(Temp, t_ind)

dudx = diff(Temp.sol(t_ind, :))./diff(Temp.x);
x = (Temp.x(1:end-1)+Temp.x(2:end))/2;
sol = (Temp.sol(t_ind, 1:end-1)+Temp.sol(t_ind, 2:end))/2;

[~, flux ,~, g0, pt] = flory_hugg_pde(x, 0, sol, dudx, Temp.a, Temp.b, ...
                Temp.e, Temp.u0, Temp.e_g0, Temp.u_g0, Temp.v, Temp.mode);
end