From f76b43f18ce467af733f33fc81977395c0a1b60d Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Mon, 8 Feb 2021 17:27:13 +0100
Subject: [PATCH] Implement skipping immediate boundary region for boundary
crossing case.
---
int_prob.m | 16 +++++++++++-----
1 file changed, 11 insertions(+), 5 deletions(-)
diff --git a/int_prob.m b/int_prob.m
index e39d8bf..32756bc 100644
--- a/int_prob.m
+++ b/int_prob.m
@@ -1,17 +1,23 @@
-function p = int_prob(l, T, x0, direc, ind_t, bp, ind_delta, T_mov)
+function p = int_prob(l, T, x0, direc, ind_t, bp, ind_delta, lb, T_mov)
% direc ... change direction of jumps, 1: left->right, -1: right->left
% 2: left->left, -2: right->right
% ind_t ... time index at which propagators are evaluated
% bp ... boundary position at t==ind_t
+% lb ... distance from boundary not to take into account
+
delta_x0 = diff(x0);
p = 0;
for i = 1:length(delta_x0)
x = (x0(i)+x0(i+1))/2;
% 1. cond.: corr. starting point? 2. cond: jumped outside of domain?
- corr_starting_point = direc*x < direc*bp;
+ corr_starting_point = direc*x < direc*bp-lb;
if abs(direc) == 1 % jump across the boundary?
- corr_end_point = direc*(x-l) > direc*(bp-T{i}.v*T{1}.t(ind_delta+1));
+ corr_end_point = direc*(x-l) > direc*(bp-T{i}.v*T{1}.t(ind_delta+1))+lb;
elseif abs(direc) == 2 % jump within the same phase?
+ if lb ~= 0
+ disp('Jump within same phase not implemented for lb!=0');
+ break;
+ end
l = abs(l); % both directions need to be considered.
if (x < bp)
if (x + l < bp); right = 1; else; right = 0; end
@@ -23,7 +29,7 @@ for i = 1:length(delta_x0)
corr_end_point = left || right; % consider, if left or right works
end
if corr_starting_point && corr_end_point
- if nargin==7
+ if nargin==8
if abs(direc) == 2
disp('Jumps within same phase not implemented for equ. .');
end
@@ -32,7 +38,7 @@ for i = 1:length(delta_x0)
% @ 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==8
+ elseif nargin==9
if abs(direc) == 1
p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_delta+1, :), x-l);
elseif abs(direc) == 2
--
GitLab