Adding first version of integrals, new formatting.

prob_laplace.m

@@ -4,32 +4,32 @@ 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 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
-                  0.2, 300, 7], t, 0.5);
+T = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3),...
+                  0.16, 0.2, 300, 7], t, 0.5);
 T.solve_tern_frap()
 %%
-T1 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
-                   0.2, 300, 7], t, 1);
+T1 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3),...
+                   0.16, 0.2, 300, 7], t, 1);
 T1.solve_tern_frap()
 %%
-T2 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
-                   0.2, 300, 7], t, 2);
+T2 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3),...
+                   0.16, 0.2, 300, 7], t, 2);
 T2.solve_tern_frap()
 %%
 tic
-T3 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
-                   0.2, 300, 7], t, 5);
+T3 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3),...
+                   0.16,0.2, 300, 7], t, 5);
 T3.solve_tern_frap()
 toc
 %%
-T4 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
-                   0.2, 300, 7], t, 10);
+T4 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3),...
+                   0.16, 0.2, 300, 7], t, 10);
 T4.solve_tern_frap()
 %% different x0
 x0 = [-0.05, -0.03, -0.01, 0.01, 0.03, 0.05];
 for i = 1:length(x0)
-T(i) = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
-                   0.2, 300, 6+x0(i)], t, 1);
+T(i) = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3),...
+                     0.16, 0.2, 300, 6+x0(i)], t, 1);
 T(i).solve_tern_frap()
 end
 %%
@@ -47,8 +47,8 @@ ga = (T(j).u0-T(j).e)/ (T(j).u0+T(j).e);
 p_out = @(D_p, D_m, ga, x0, x, t) 1./(2*sqrt(D_p*pi*t))*...
             (exp(-(x+x0).^2./(4*D_p*t))*(ga*sqrt(D_p)-sqrt(D_m))./...
             (ga*sqrt(D_p)+sqrt(D_m))+exp(-(x-x0).^2./(4*D_p*t)))*2;
-p_in = @(D_p, D_m, ga, x0, x, t) 1./(sqrt(pi*t)*(sqrt(D_m)+ga*sqrt(D_p)))*...
-            exp(-(x-x0*sqrt(D_m/D_p)).^2/(4*D_m*t))*2;
+p_in = @(D_p, D_m, ga, x0, x, t) 1./(sqrt(pi*t)*(sqrt(D_m)+...
+            ga*sqrt(D_p)))*exp(-(x-x0*sqrt(D_m/D_p)).^2/(4*D_m*t))*2;
 x_left = linspace(-4, 0, 1000);
 x_right = linspace(0, 4, 1000);
 %% Plot with full ternary model
@@ -79,4 +79,44 @@ for i = 1:4
 pks_min = findpeaks(-T4.sol(i, :));
 pks_max = findpeaks(T4.sol(i, :));
 pks_max(1)/pks_min(1)
+end
+
+%% Solve integrals
+
+% different x0
+clear T;
+x0 = 0.01:0.01:10;
+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];
+% disp([-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
+%                    0.2, 300, 6+x0(i)]);
+T{i} = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), ...
+                     0.16, 0.2, 300, 6+x0(i)], t, 1);
+T{i}.solve_tern_frap()
+end
+toc
+%%
+ls = 0.01:0.01:2;
+for i = 1:length(ls)
+q(i) = int_prob(ls(i), T, x0+6);
+end
+%%
+figure; hold on;
+plot(ls, p);
+plot(ls, q);
+%% distribution for x0 can be taken from phi_tot (steady state)
+function p = int_prob(l, T, x0)
+delta_x0 = diff(x0);
+p = 0;
+for i = 1:length(delta_x0)
+    x = (x0(i)+x0(i+1))/2;
+    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 = p + delta_x0(i)*...
+        T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0)  * p2;
+    if x0(i)-l < 0; break; end
+end
 end
\ No newline at end of file