Simplified probability integration. Goes to zero. length scales wrong.

prob_laplace.m

@@ -82,38 +82,53 @@ pks_max(1)/pks_min(1)
 end
 
 %% Solve integrals
-
-% different x0
-clear T;
-x0 = 0.0:0.01:2.02;
+a = 5;
+nu = 10^-6;
+chi = 7/3;
+b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
+e = @(chi) sqrt(3/8*(chi-2));
+t = [0, 0.05, 0.1];
+T1 = Ternary_model(0, 'Gauss', [-5, b(chi, nu), 0.5, e(chi), ...
+                     0.0, 0.0, 10, a], t, 3);
+x0 = T1.x(find(T1.x>a, 1):end);
+x0 = x0(1:10:end);
+% x0 = 0.0:0.001:0.1;
+%%
+T = {};
 tic
 parfor i = 1:length(x0)
 b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
 e = @(chi) sqrt(3/8*(chi-2));
-t = [0, 0.1, 1, 9.9];
-T{i} = Ternary_model(0, 'Gauss', [-5, b(7.7/3, 10^-6), 0.5, e(7.7/3), ...
-                     0.16, 0.2, 10, 5+x0(i)], t, 3);
+T{i} = Ternary_model(0, 'Gauss', [-a, b(chi, nu), 0.5, e(chi), ...
+                     0.0, 0.0, 10, x0(i)], t, 3);
 T{i}.solve_tern_frap()
 end
 toc
 %%
-ls = 0.0:0.0005:1;
-for i = 1:length(ls)
-q(i) = int_prob(ls(i), T, x0+5);
+ls = 0.001:0.01:1;
+q = nan(1, length(ls));
+parfor i = 1:length(ls)
+q(i) = int_prob_simple(ls(i), T, x0);
 end
-%%
-Ternary_model(0, 'Gauss', [-5, b(7.7/3, 10^-6), 0.5, e(7.7/3), ...
-                     0.16, 0.2, 10, 6+x0(i)], t, 3)
 %% Do we need to look at the left side as well?
 figure; hold on;
 % plot(ls, p);
-plot(ls, q);
+plot(ls, N*q);
+%%
+int_prob(0.01, T, x0+a)
 %%
-int_prob(0.0, T, x0+5)
+int_prob_simple(0.01, T, x0)
 %%
 for i = 1:7:length(T)
     T{i}.plot_sim('plot', 1, 'red')
 end
+%% Normalization factor
+tic
+N = normalization(T, x0+a);
+toc
+%% Write to file
+csvwrite('jump_length.csv', ls)
+csvwrite('prob.csv', q);
 %% distribution for x0 can be taken from phi_tot (steady state)
 function p = int_prob(l, T, x0)
 delta_x0 = diff(x0);
@@ -121,9 +136,49 @@ p = 0;
 for i = 1:length(delta_x0)
     x = (x0(i)+x0(i+1))/2;
     if x-l < 0; break; end
-    p_i = @(j, x) interp1(T{j}.x, T{j}.sol(2, :), x-l);
-    p2 = (p_i(i, x)+p_i(i+1, x))/2;
+    p_i = @(j) interp1(T{j}.x, T{j}.sol(2, :), 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;
 end
 end
+
+function p = int_prob_simple(l, T, x0)
+p = 0;
+p_i = @(j, x0) interp1(T{j}.x, T{j}.sol(3, :), x0-l);
+for i = 1:length(x0)
+    if x0(i)-l > 5; break; end % Is this right?
+    % Calculate left bin boundary
+    if i == 1
+        left_bound = -T{1}.a; 
+    else
+        left_bound = (x0(i)-x0(i-1))/2;
+    end
+    % calculate right bin boundary
+    if i == length(x0)
+        right_bound = T{1}.system_size;
+    else
+        right_bound = (x0(i+1)+x0(i))/2;
+    end
+    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)*...
+            p_i(i, x0(i));
+end
+end
+
+function p = normalization(T, x0)
+delta_x0 = diff(x0);
+p = 0;
+parfor i = 1:length(delta_x0)
+    x = (x0(i)+x0(i+1))/2;
+    p_x0i = T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0);
+    for j = 1:length(delta_x0)
+        xj = (x0(j)+x0(j+1))/2;
+        p_j = @(j) interp1(T{i}.x, T{i}.sol(2, :), xj);
+        pj = (p_j(i)+p_j(i+1))/2;
+        p = p + delta_x0(i)*delta_x0(j)*p_x0i*pj;
+    end
+end
+end
\ No newline at end of file