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

ss = T{1}.system_size;
delta_x0 = diff(x0);
p = 0;

if lb ~= 0 && abs(direc)==2
    disp('Jump within same phase not implemented for lb!=0');
    p = nan;
    return;
end

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-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))+lb;
    elseif direc == 2 % jump left -> left
        if l < 0 && x-l < bp; corr_end_point = 1; else; corr_end_point=0; end
        if l > 0 % jump to left
            corr_end_point = 1;
            if x-l < 0; x = l-x; end % reflecting boundary if jump out left
        end
    elseif direc == -2 % jump right -> right
        if l > 0 && x-l > bp; corr_end_point = 1; else; corr_end_point=0; end
        if l < 0 % jump to right
            corr_end_point = 1;
            % Reflecting boundary if jumping out on right side of system
            if x-l > ss; x = 2*ss+l-x; end
        end
    end
    if corr_starting_point && corr_end_point
        if nargin==8
            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==9
            if abs(direc) == 1
                x_ev = x-l;
            elseif direc == 2 % left -> left
                if l > 0 && x-l < 0; x_ev = l-x; else; x_ev = x-l; end
            elseif direc == -2 % right -> right
                if l < 0 && x+l > ss; x_ev = 2*ss+l-x; else; x_ev = x-l; end
            end
            p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_delta+1, :), x_ev);
            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