%% 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, '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
%% %%%%%%%%%%%%%%%%% 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', 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 = 4;
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);
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_lb01.csv', ls)
csvwrite('prob_7_7_lb01.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', 0.5, 10, 7, 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/';
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', 0.5, 10, 7, 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/';
csvwrite([pa, 'jump_length_7_7_lb002_bp', num2str(bp(j)), '.csv'], ls)
csvwrite([pa, 'prob_7_7_lb002_bp', num2str(bp(j)), '.csv'], p);
end
%% Mean cross jump length chi 7/3
run_jump_lengths(-5, 1, 7/3, 1, 4, '73_1_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, 1, 7/3, 1, 4, '73_1_10', 10^-6, 4, 0.5);
run_jump_lengths(-5, 0, 7/3, 1, 4, '73_6_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, 0, 7/3, 1, 4, '73_6_10', 10^-6, 4, 0.5);
run_jump_lengths(-5, -0.08410041, 7/3, 1, 4, '73_10_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, -0.08410041, 7/3, 1, 4, '73_10_10', 10^-6, 4, 0.5);
run_jump_lengths(-5, -0.12314694, 7/3, 1, 4, '73_15_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, -0.12314694, 7/3, 1, 4, '73_15_10', 10^-6, 4, 0.5);
run_jump_lengths(-5, -0.14264143, 7/3, 1, 4, '73_20_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, -0.14264143, 7/3, 1, 4, '73_20_10', 10^-6, 4, 0.5);
run_jump_lengths(-5, -0.15432893, 7/3, 1, 4, '73_25_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, -0.15432893, 7/3, 1, 4, '73_25_10', 10^-6, 4, 0.5);
run_jump_lengths(-5, -0.16211677, 7/3, 1, 4, '73_30_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, -0.16211677, 7/3, 1, 4, '73_30_10', 10^-6, 4, 0.5);
%% chi 7.7/3
run_jump_lengths(-5, 1, 7.7/3, 1, 4, '773_1_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, 1, 7.7/3, 1, 4, '773_1_10', 10^-6, 4, 0.5);
run_jump_lengths(-5, 0.136322219, 7.7/3, 1, 4, '773_6_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, 0.136322219, 7.7/3, 1, 4, '773_6_10', 10^-6, 4, 0.5);
run_jump_lengths(-5, 0.0618145816, 7.7/3, 1, 4, '773_10_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, 0.0618145816, 7.7/3, 1, 4, '773_10_10', 10^-6, 4, 0.5);
run_jump_lengths(-5, 0.027120457, 7.7/3, 1, 4, '773_15_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, 0.027120457, 7.7/3, 1, 4, '773_15_10', 10^-6, 4, 0.5);
run_jump_lengths(-5, 0.00977482501, 7.7/3, 1, 4, '773_20_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, 0.00977482501, 7.7/3, 1, 4, '773_20_10', 10^-6, 4, 0.5);
run_jump_lengths(-5, 0, 7.7/3, 1, 4, '773_25_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, 0, 7.7/3, 1, 4, '773_25_10', 10^-6, 4, 0.5);
run_jump_lengths(-5, -0.00756985349, 7.7/3, 1, 4, '773_30_01', 10^-6, 3, 0.5);
run_jump_lengths(-5, -0.00756985349, 7.7/3, 1, 4, '773_30_10', 10^-6, 4, 0.5);
% output is saved manually from displayed values to data/mean_cross_jl_ML.csv
%% %%%%%%%%%%%%%%%%%%% 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, p_i_i] = frap_distro(t, params, 2, bp, 1);
[~, p_o_o] = frap_distro(t, params, -2, bp, 1);
[~, p_i_o] = frap_distro(t, params, 1, bp, 1);
[~, 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;
%%
figure; hold on;
plot(ls, sum(p_i_i, 2));
plot(ls, sum(p_o_o, 2));
plot(ls, sum(p_i_o, 2));
plot(ls, sum(p_o_i, 2));
save('/data/biophys/hubatsch/MatlabWorkspaces/prob_laplace_190121.mat');
%% Save time points for comparison with python
csvwrite('lb_in_out.csv', ls)
csvwrite('prob_in_out.csv', p);
%%
KL = zeros(1, 100);
for i = 1:68
i_t = 1+i;
j_t = 2+i;
integrand = @(p, q) abs(p.*log(p./q));
integrate = @(x, y) nansum(diff(x)'.*(y(2:end)+y(1:end-1))/2);
n_i = @(p1, p2, p3, p4) abs(integrate(ls, p1) + integrate(ls, p2)+...
integrate(ls, p3) + integrate(ls, p4));
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));
KL(i) = integrate(ls, integrand(p_i_i(:, i_t)/N_i, p_i_i(:, j_t)/N_j))+...
integrate(ls, integrand(p_o_o(:, i_t)/N_i, p_o_o(:, j_t)/N_j))+...
integrate(ls, integrand(p_i_o(:, i_t)/N_i, p_o_i(:, j_t)/N_j))+...
integrate(ls, integrand(p_o_i(:, i_t)/N_i, p_i_o(:, j_t)/N_j));
end
plot(t(2:101), abs(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.1908/0.1885*0.61*abs(KL)/0.02);
%% %%%%%%%%%%%%%%%%%%% 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);
n_T=70;
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