%% 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.01, 0.05, 50];
%% 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)
    tic
    b = 2e-4;
    u0 = 0.5;
    P=1;
    D_i = 2;
    D_o = 20;
    a = -1;
    [b, u0, 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, 2000);
%     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()
    toc
end

%% Test partitioning == 1 compared to natural no mobility case.
% Tg1: g0 = 1 everywhere. 

Tpt1 = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-18), 0.5, e(7/3),...
                                0, 1, 300, 7, 0, 'Constituent'},...
                     t, 1);
Tpt1.solve_tern_frap();
%%
% Rescale mobility, such that phi_tot=1 everywhere achieves same flux
Gi = (1-Tpt1.phi_t(1))/(1-Tpt1.phi_t(1));
Go = Tpt1.phi_t(1)/Tpt1.phi_t(end)*(1-Tpt1.phi_t(end));
%%
Tga1 = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-18), Tpt1.phi_t(1), 0,...
                                 (Go-Gi), Gi, 300, 7, 0, 'Constituent'},...
                                 t, 1);
Tga1.solve_tern_frap();
%%
Tpt1.plot_sim('plot', 1, 'g')
%%
Tga1.plot_sim('plot', 1, 'k')
%%
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, 'Constituent'}, 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
%% Moving boundary
t = linspace(0, 10, 9000);
tic
T_mov = Ternary_model(0, 'phi_tot', {-10, b(7/3, 10^-6), 0.5, e(7/3),...
                   0, 1, 10, 7, 0.5, 'Client'}, t, 0);
T_mov.solve_tern_frap();
toc
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;
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;
g0 = 0.5;
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);
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])
%% 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, 'Constituent'}, 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));
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


%% Solve integrals for jump length distribution @ steady state.
% Set parameters
params = {-5, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0, 1, 10, 7, 0, 'Constituent'};
t = [0, 0.05, 0.1];
x0 = 5.0:0.002:7.01;
%% Run simulations with '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.
ls = 0.0005:0.005:2;
p = nan(1, length(ls));
tic
parfor i = 1:length(ls)
p(i) = int_prob(ls(i), T, x0, 0);
end
toc
%% Test integrals
int_prob(0.3, T, x0, 0.1)
int_prob_simple(0.01, T, x0)
%% Normalization factor
N = normalization(T, x0, 0);
%% Do we need to look at the left side as well?
figure; hold on;
plot(ls, p/N);
%% 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, T_mov, ind_t, ind_delta)
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
        if nargin==4
            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, 0)*p2;
        elseif nargin==7
            p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_t+ind_delta, :), x-l);
            p2 = (p_i(i)+p_i(i+1))/2;
            p_sol = @(j) interp1(T_mov.x, T_mov.sol(ind_t, :), x);
            p_sol2 = (p_sol(i)+p_sol(i+1))/2;
            p = p + delta_x0(i)*p_sol2*p2;
        end
    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, 0)*...
            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, 0);
        p_i(i) = delta_x0(i)*sum(sol.*delta_x)*p_x0i;
    end
end
p = sum(p_i);
end