%% 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.05, 1, 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^-9), 0.5, e(7/3),...
                                0, 1, 300, 7, 0, 'Constituent'},...
                     t, 0.5);
% Tpt1.x = Tga1.x;
Tpt1.solve_tern_frap();
%%
% 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;
%%
Tga1 = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-9), Tpt1.phi_t(1), 0,...
                                 0, Gi, 300, 7, 0, 'Constituent'},...
                                 t, 0.5);
% Tga1.x = Tpt1.x;
Tga1.solve_tern_frap();
%%
Tpt1.plot_sim('plot', 1, 'r', 1.72)
Tga1.plot_sim('plot', 1, 'k', 1)

%% Calculate flux for Tpt1 at t=0
fac = 1;
[dudx, x, sol, flux, g0, pt] = calc_flux(Tpt1, 1);
[dudx1, x1, sol1, flux1, g01, pt1] = calc_flux(Tga1, 1);
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, '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 TIME COURSE %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 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'}, 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


%% %%%%%%%%%%%%%%%%% STEADY STATE JUMP LENGTH DISTRIBUTION %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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, 1];
direction = -1;
x0 = sort(5-direction*(0:0.002:3.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.
ls = 0.0005:0.005:3;
p = nan(1, length(ls));
tic
parfor i = 1:length(ls)
    p(i) = int_prob(ls(i), T, x0, direction, 4);
end
toc
%% Test integrals
int_prob(0.3, T, x0, direction, 4)
%% Normalization factor
N = normalization(T, x0, 0, 4);
%% Do we need to look at the left side as well?
figure; hold on;
plot(ls, p/N);
%% Write to file
cd /Users/hubatsch/Nextcloud/Langevin_vs_MeanField/Data_Figs_FokkPla/jump_length/
csvwrite('jump_length_7_7_lb01.csv', ls)
csvwrite('prob_7_7_lb01.csv', p);


%% %%%%%%%%%%%%%%%%%%%%%% FRAP JUMP LENGTH DISTRIBUTION %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
params = {-5, b(7/3, 10^-6), 0.5, e(7/3), 0, 1, 10, 7, 0, 'Constituent'};
t = linspace(0, 5, 51);
direction = 1; % 1: jump from left to right, -1: jump from right to left.
x0 = sort(5-direction*(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 = -direction*(0.001:0.04:4);
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, direction, T_mov, j, 5);
    end
    toc
end
% save prob_laplace_X_7_3_FRAP_in_out
%% Save time points for comparison with python
csvwrite('lb_in_out.csv', ls)
csvwrite('prob_in_out.csv', p);

%% %%%%%%%%%%%%%%%%%%% INTEGRATION FUNCTION DEFINITIONS %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = int_prob(l, T, x0, direction, ind_t, T_mov, ind_delta)
% 'direction' change direction of jumps, 1: left->right, -1: right->left
delta_x0 = diff(x0);
p = 0;
for i = 1:length(delta_x0)
    x = (x0(i)+x0(i+1))/2;
    if (direction*x < direction*5) && (direction*(x-l) > direction*5)
        if nargin==5
            p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_t, :), x-l);
            p2 = (p_i(i)+p_i(i+1))/2;
            % @ SS distribution for x0 can be taken from phi_tot
            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_delta+1, :), 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 = T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0, 0);
            p_sol2 = (p_sol(i)+p_sol(i+1))/2;
            p = p + delta_x0(i)*p_sol2*p2;
        end
    end
end
end

function p = normalization(T, x0, lb, ind_t)
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(ind_t, 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

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