Skip to content
Snippets Groups Projects
Commit 5d47b16e authored by Lars Hubatsch's avatar Lars Hubatsch
Browse files

Jump distribution work bidirectionally.

parent 08783414
No related branches found
No related tags found
No related merge requests found
......@@ -27,7 +27,7 @@ elseif strcmp(mode, 'plot')
figure(20);
cla;
hold on;
xlim([-inf, 1.0]); ylim([-inf, max(T.sol(:))]);
xlim([-inf, 6.0]); ylim([-inf, max(T.sol(:))]);
ax = gca;
ax.FontSize = 12;
xlabel('x [\mum]'); ylabel('volume fraction');
......
......@@ -205,8 +205,9 @@ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
direction = -1;
x0 = sort(5-direction*(0:0.002:2.01));
%% Run simulations for 'delta' IC across outside
T = {};
parfor i = 1:length(x0)
tic
......@@ -221,11 +222,11 @@ 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);
p(i) = int_prob(ls(i), T, x0, direction);
end
toc
%% Test integrals
int_prob(0.3, T, x0, 0.1)
int_prob(0.3, T, x0, direction)
int_prob_simple(0.01, T, x0)
%% Normalization factor
N = normalization(T, x0, 0);
......@@ -242,8 +243,9 @@ csvwrite('prob_7_7_lb01.csv', p);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
params = {-5, b(7/3, 10^-6), 0.5, e(7/3), 0, 1, 10, 7, 0, 'Constituent'};
t = linspace(0, 5, 51);
x0 = 5.0:0.002:9.01;
%% Run simulations with 'delta' IC across outside
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
......@@ -255,32 +257,30 @@ 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 = 0.001:0.04:4;
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, 0, T_mov, j, 5);
p(i, j) = int_prob(ls(i), F, x0, direction, T_mov, j, 5);
end
toc
end
save prob_laplace_X_7_3_FRAP
%%
int_prob(0.05, F, x0, 0, T_mov, 1, 30)
% save prob_laplace_X_7_3_FRAP_in_out
%% Save time points for comparison with python
csvwrite('lb.csv', ls)
csvwrite('prob.csv', p);
csvwrite('lb_in_out.csv', ls)
csvwrite('prob_in_out.csv', p);
%% %%%%%%%%%%%%%%%%%%% INTEGRATION FUNCTION DEFINITIONS %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = int_prob(l, T, x0, lb, T_mov, ind_t, ind_delta)
function p = int_prob(l, T, x0, direction, T_mov, ind_t, 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 x-5>lb
if x-l > 5-lb; break; end
if (direction*x < direction*5) && (direction*(x-l) > direction*5)
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;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment