Change signature of int_prob. SS jump lengths work in both directions.

prob_laplace.m

@@ -189,10 +189,11 @@ 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'};
+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.002:3.01));
+direc = 1;
+x0 = sort(5-direc*(0:0.001:3.01));
 %% Run simulations for 'delta' IC across outside
 T = {};
 parfor i = 1:length(x0)
@@ -204,15 +205,17 @@ parfor i = 1:length(x0)
 end
 % save prob_laplace_X_7_7_short_2
 %% Calculate probabilities for each jump length in ls.
-ls = 0.0005:0.005:3;
+ls = -direc*(0.000:0.001:3);
 p = nan(1, length(ls));
 tic
 parfor i = 1:length(ls)
-    p(i) = int_prob(ls(i), T, x0, direc, 4, 5);
+    tic
+    p(i) = int_prob(ls(i), T, x0, direc, 3, 5, 0);
+    toc
 end
 toc
 %% Normalization factor
-N = normalization(T, x0, 0, 4);
+N = normalization(T, x0, 0, 3);
 %% Do we need to look at the left side as well?
 figure; hold on;
 plot(ls, p/N);
@@ -293,7 +296,7 @@ for j = 1:n_T
     F1 = F(:, j);
     bp = bound(j);
     for i = 1:length(ls)
-        p(i, j) = int_prob(ls(i), F1, x0(j, :), direc, j, bp, T_mov, 10);
+        p(i, j) = int_prob(ls(i), F1, x0(j, :), direc, j, bp, 10, T_mov);
     end
     toc
 end
@@ -377,7 +380,7 @@ 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, T_mov, 5);
+        p(i, j) = int_prob(ls(i), F, x0, direc, j, bp, 5, T_mov);
     end
     toc
 end
@@ -388,7 +391,7 @@ csvwrite('prob_in_out.csv', p);
 
 %% %%%%%%%%%%%%%%%%%%% INTEGRATION FUNCTION DEFINITIONS %%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-function p = int_prob(l, T, x0, direc, ind_t, bp, T_mov, ind_delta)
+function p = int_prob(l, T, x0, direc, ind_t, bp, ind_delta, T_mov)
 % direc ... change direction of jumps, 1: left->right, -1: right->left
 % ind_t     ... time index at which propagators are evaluated
 % bp        ... boundary position at t==ind_t
@@ -400,7 +403,7 @@ for i = 1:length(delta_x0)
     corr_starting_point = direc*x < direc*bp;
     corr_end_point = direc*(x-l) > direc*(bp-T{i}.v*T{1}.t(ind_delta+1));
     if corr_starting_point && corr_end_point
-        if nargin==6
+        if nargin==7
             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