Outsourcing jump length calculation to separate function.

prob_laplace.m

@@ -185,25 +185,40 @@ for i = 1:100%length(T_mov.t)
 %     print([num2str(i),'.png'],'-dpng')
     shg; pause();
 end
+%% 
+run_jump_lengths(-5, -0.15432893, 7/3, 1, 4, '73_0154_10', 10^-6, 4, 0.5);
+%% chi 7/3
+run_jump_lengths(-5, 1, 7/3, 1, 4, '73_1_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, 1, 7/3, 1, 4, '73_1_10', 10^-6, 4, 0.5);
+run_jump_lengths(-5, 0, 7/3, 1, 4, '73_6_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, 0, 7/3, 1, 4, '73_6_10', 10^-6, 4, 0.5);
+run_jump_lengths(-5, -0.08410041, 7/3, 1, 4, '73_10_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, -0.08410041, 7/3, 1, 4, '73_10_10', 10^-6, 4, 0.5);
+run_jump_lengths(-5, -0.12314694, 7/3, 1, 4, '73_15_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, -0.12314694, 7/3, 1, 4, '73_15_10', 10^-6, 4, 0.5);
+run_jump_lengths(-5, -0.14264143, 7/3, 1, 4, '73_20_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, -0.14264143, 7/3, 1, 4, '73_20_10', 10^-6, 4, 0.5);
+run_jump_lengths(-5, -0.15432893, 7/3, 1, 4, '73_25_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, -0.15432893, 7/3, 1, 4, '73_25_10', 10^-6, 4, 0.5);
+run_jump_lengths(-5, -0.16211677, 7/3, 1, 4, '73_30_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, -0.16211677, 7/3, 1, 4, '73_30_10', 10^-6, 4, 0.5);
 
-%% Different ratios of Dm/Dp
-chi = 7/3; nu = 10^-6; u0 = 0.5; a = -5;
-b1 = b(chi, nu); P = (u0+e(chi))/(u0-e(chi));
-params = {a, b1, u0, e(chi), -0.15432893, 0, 10, 7, 0, 'Const_mob', 0};
-%% Stefano's way of having mobility with square. works, 0.5% diff
-chi = 7/3; nu = 10^-6; u0 = 0.5; a = -5;
-b1 = b(chi, nu); P = (u0+e(chi))/(u0-e(chi));
-params = {a, b1, u0, e(chi), -0.12314694, 0, 10, 7, 0, 'Const_mob', 0};
-%%
-D_i = 0.05672749; D_o = 0.85091231; chi = 7/3; nu = 10^-6; u0 = 0.5; a = -5;
-b1 = b(chi, nu); P = (u0+e(chi))/(u0-e(chi));
-[~, e_g, u_g] = calc_tanh_params(P, D_i, D_o, a, b1, u0);
-%%
-D_i = 0.03; D_o = 0.83; chi = 7/3; nu = 10^-6; u0 = 0.5; a = -5;
-b1 = b(chi, nu); P = (u0+e(chi))/(u0-e(chi));
-[~, e_g, u_g] = calc_tanh_params(P, D_i, D_o, a, b1, u0);
-%%
-params = {a, b1, u0, e(chi), e_g, u_g, 10, 7, 0, 'Constituent', 0};
+replace(num2str(round(chi,2)),'.', '')
+%% chi 7.7/3
+run_jump_lengths(-5, 1, 7.7/3, 1, 4, '773_1_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, 1, 7.7/3, 1, 4, '773_1_10', 10^-6, 4, 0.5);
+run_jump_lengths(-5, 0.136322219, 7.7/3, 1, 4, '773_6_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, 0.136322219, 7.7/3, 1, 4, '773_6_10', 10^-6, 4, 0.5);
+run_jump_lengths(-5, 0.0618145816, 7.7/3, 1, 4, '773_10_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, 0.0618145816, 7.7/3, 1, 4, '773_10_10', 10^-6, 4, 0.5);
+run_jump_lengths(-5, 0.027120457, 7.7/3, 1, 4, '773_15_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, 0.027120457, 7.7/3, 1, 4, '773_15_10', 10^-6, 4, 0.5);
+run_jump_lengths(-5, 0.00977482501, 7.7/3, 1, 4, '773_20_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, 0.00977482501, 7.7/3, 1, 4, '773_20_10', 10^-6, 4, 0.5);
+run_jump_lengths(-5, 0, 7.7/3, 1, 4, '773_25_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, 0, 7.7/3, 1, 4, '773_25_10', 10^-6, 4, 0.5);
+run_jump_lengths(-5, -0.00756985349, 7.7/3, 1, 4, '773_30_01', 10^-6, 3, 0.5);
+run_jump_lengths(-5, -0.00756985349, 7.7/3, 1, 4, '773_30_10', 10^-6, 4, 0.5);
 %% %%%%%%%%%%%%%%%%% STEADY STATE JUMP LENGTH DISTRIBUTION %%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 params = {-5, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0, 1, 10, 7, 0,...
@@ -235,7 +250,8 @@ parfor i = 1:length(ls)
 end
 toc
 %% Normalization factor
-N = normalization(T, x0, 0, t_ind, direc, -params{1})
+% N = normalization(T, x0, 0, t_ind, direc, -params{1})
+N=sum(p)/length(p)*l_max
 % N = sum(p)*0.001;
 m = sum(ls.*p/N)/length(p)*l_max;
 % should sum to one

run_jump_lengths.m

+function run_jump_lengths(a, ad, chi, direc, l_max, name, nu, t_ind, u0)
+% Calculate jump length distribution for given parameter set.
+% ad ... prefactor of quadratic term in Stefano's mobility ansatz
+b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
+e = @(chi) sqrt(3/8*(chi-2));
+params = {a, b(chi, nu), u0, e(chi), ad, 0, 10, 7, 0, 'Const_mob', 0};
+t = [0, 0.05, 0.1, 1];
+x0 = sort(5-direc*(0:0.001:4.01));
+%% Run simulations for 'delta' IC across outside
+T = {};
+parfor i = 1:length(x0)
+    tic
+    T{i} = Ternary_model(0, 'Gauss', params, t, 0.2);
+    T{i}.x0 = x0(i);
+    T{i}.solve_tern_frap();
+    toc
+end
+%% Calculate probabilities for each jump length in ls.
+ls = -direc*(0.000:0.001:l_max);
+p = nan(1, length(ls));
+tic
+parfor i = 1:length(ls)
+    tic
+    p(i) = int_prob(ls(i), T, x0, direc, t_ind, 5, 0);
+    toc
+end
+toc
+%% Normalization factor
+% N = normalization(T, x0, 0, t_ind, direc, -params{1});
+% should sum to one
+N = sum(p)/length(p)*l_max;
+m = sum(ls.*p/N)/length(p)*l_max;
+save(name);