From f630c23fd90456149c6fd565163680a3bf796366 Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Tue, 19 Jan 2021 14:30:00 +0100
Subject: [PATCH] out -> out and in -> in corrected.

---
 int_prob.m | 26 ++++++++++++++++++++++++--
 1 file changed, 24 insertions(+), 2 deletions(-)

diff --git a/int_prob.m b/int_prob.m
index 60fc26c..e39d8bf 100644
--- a/int_prob.m
+++ b/int_prob.m
@@ -12,17 +12,39 @@ for i = 1:length(delta_x0)
     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));
     elseif abs(direc) == 2 % jump within the same phase?
-        corr_end_point = direc*(x-l) < direc*(bp-T{i}.v*T{1}.t(ind_delta+1));
+        l = abs(l); % both directions need to be considered.
+        if (x < bp)
+            if (x + l < bp);    right = 1;  else; right = 0;    end
+            if (x - l > 0);   left = 1;   else; left = 0;     end
+        else
+            if (x + l < T{1}.system_size);  right =1;   else; right = 0;    end
+            if (x -l > bp);     left = 1;   else; left = 0;     end
+        end
+        corr_end_point = left || right; % consider, if left or right works
     end
     if corr_starting_point && corr_end_point
         if nargin==7
+            if abs(direc) == 2
+                disp('Jumps within same phase not implemented for equ. .');
+            end
             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==8
-            p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_delta+1, :), x-l);
+            if abs(direc) == 1
+                p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_delta+1, :), x-l);
+            elseif abs(direc) == 2
+                x_l = x-l; % left position
+                x_r = x+l; % right position
+                if x-l < 0; x_temp = l-x; end % -> jump reflected if necess.
+                if x+l > T{1}.system_size; x_temp = 2*T{1}.system_size-x-l; end
+                p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_delta+1, :),...
+                                   x_l) * left + ...
+                           interp1(T{j}.x, T{j}.sol(ind_delta+1, :),...
+                                   x_r) * right;
+            end
             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);
-- 
GitLab