Preparing prob_laplace to add FRAP integrals. Tested against prob_7_7_lb01 (Apr.).

prob_laplace.m

@@ -192,66 +192,63 @@ pks_max = findpeaks(T_prec(5).sol(i, :));
 pks_max(1)/pks_min(1)
 end
 
-%% Solve integrals
-a = 5;
-nu = 10^-6;
-chi = 7.7/3;
-b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
-e = @(chi) sqrt(3/8*(chi-2));
+
+%% Solve integrals for jump length distribution @ steady state.
+% Set parameters
+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];
-% T1 = Ternary_model(0, 'Gauss', [-5, b(chi, nu), 0.5, e(chi), ...
-%                      0, 1, 10, a], t, 0.2);
 x0 = 5.0:0.002:7.01;
-%%
+%% Run simulations with 'delta' IC across outside
 T = {};
 parfor i = 1:length(x0)
 tic
-b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
-e = @(chi) sqrt(3/8*(chi-2));
-T{i} = Ternary_model(0, 'Gauss', [-a, b(chi, nu), 0.5, e(chi), ...
-                     0, 1, 10, x0(i)], t, 0.2);
+T{i} = Ternary_model(0, 'Gauss', params, t, 0.2);
+T{i}.x0 = x0(i);
 T{i}.solve_tern_frap();
 toc
 end
 % save prob_laplace_X_7_7_short_2
-%%
+%% Calculate probabilities for each jump length in ls.
 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.1);
+p(i) = int_prob(ls(i), T, x0, 0);
 end
-%% Do we need to look at the left side as well?
-figure; hold on;
-plot(ls, p/N);
-% plot(lx, xt);
-% plot(ls, q/50000);
-%%
-for i = 1:2:length(T)
-    T{i}.plot_sim('plot', 3, 'red')
-end
-%%
+toc
+%% Test integrals
 int_prob(0.3, T, x0, 0.1)
-%%
 int_prob_simple(0.01, T, x0)
 %% Normalization factor
-tic
-N = normalization(T, x0,0.1);
-toc
+N = normalization(T, x0, 0);
+%% Do we need to look at the left side as well?
+figure; hold on;
+plot(ls, p/N);
 %% Write to file
 csvwrite('jump_length_7_7_lb01.csv', ls)
 csvwrite('prob_7_7_lb01.csv', p/N);
-%% distribution for x0 can be taken from phi_tot (steady state)
-function p = int_prob(l, T, x0, lb)
+
+
+%% Distribution for x0 can be taken from phi_tot (steady state)
+function p = int_prob(l, T, x0, lb, T_mov, ind_t, ind_delta)
 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
-        p_i = @(j) interp1(T{j}.x, T{j}.sol(3, :), x-l);
-        p2 = (p_i(i)+p_i(i+1))/2;
-        p = p + delta_x0(i)*...
-                T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0)*p2;
+        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;
+            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==7
+            p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_t+ind_delta, :), x-l);
+            p2 = (p_i(i)+p_i(i+1))/2;
+            p_sol = @(j) interp1(T_mov.x, T_mov.sol(ind_t, :), x);
+            p_sol2 = (p_sol(i)+p_sol(i+1))/2;
+            p = p + delta_x0(i)*p_sol2*p2;
+        end
     end
 end
 end
@@ -276,7 +273,7 @@ for i = 1:length(x0)
     bin_width = right_bound-left_bound;
     
     p = p + bin_width*...
-            T{1}.phi_tot(x0(i), T{1}.a, T{1}.b, T{1}.e, T{1}.u0)*...
+            T{1}.phi_tot(x0(i), T{1}.a, T{1}.b, T{1}.e, T{1}.u0, 0)*...
             p_i(i, x0(i));
 end
 end
@@ -291,7 +288,7 @@ for i = 1:length(delta_x0)
     sol = (sol(1:end-1)+sol(2:end))/2;
     x = (x0(i)+x0(i+1))/2;
     if x > 5+lb
-        p_x0i = T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0);
+        p_x0i = T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0, 0);
         p_i(i) = delta_x0(i)*sum(sol.*delta_x)*p_x0i;
     end
 end