WIP on integrating probabilites for jump length. This roughly works, but not perfect.

@Ternary_model/create_mesh.m

@@ -10,20 +10,22 @@ while x<-2*T.a
     phi_nm1 = T.phi_tot(x(end), T.a, T.b, T.e, T.u0);
     phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0);
     if x(end) <= -T.a
-        while phi_nm1-phi_n > 0.001/T.precision
+        while phi_nm1-phi_n > 0.0002/T.precision
             x_temp = x_temp/2;
             phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0);
         end
     elseif x(end) > -T.a
-        while (phi_nm1-phi_n < 0.0005/T.precision) && (x_temp < 0.01/T.precision)
+        while (phi_nm1-phi_n < 0.0001/T.precision) && ...
+                (x_temp < 0.01/T.precision)
             x_temp = 2*x_temp;
             phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0);
         end
     end
     x = [x, x(end)+x_temp];
 end
-x_middle = linspace(-2*T.a+x(end)-x(end-1), -10*T.a, 100*T.precision);
-x_right = linspace(-10*T.a+0.1, T.system_size, 100*T.precision);
+x_middle = linspace(-2*T.a+x(end)-x(end-1), -10*T.a, 200*T.precision);
+x_right = linspace(-10*T.a+0.1, T.system_size, 200*T.precision);
 T.x = [x, x_middle, x_right];
+T.x = T.x(1:find(T.x-T.system_size>=0, 1)-1);
 end
 

@Ternary_model/plot_sim.m

@@ -29,9 +29,10 @@ elseif strcmp(mode, 'plot')
     ax = gca;
     ax.FontSize = 26;
     xlabel('x [\mum]'); ylabel('volume fraction'); 
-    plot(T.x, T.sol(1:nth:si(1), :), 'LineWidth', 1, 'Color', color);%[0.5294, 0.8, 0.98]);%
-    plot(T.x, Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0), 'LineWidth', 1,...
-    'LineStyle', '--', 'Color', [0.97, 0.55, 0.03]);
+%     plot(T.x, T.sol(1:nth:si(1), :), 'LineWidth', 1, 'Color', color);
+    plot(T.x, T.sol(nth, :), 'LineWidth', 1, 'Color', color);
+    plot(T.x, Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0),...
+        'LineWidth', 1, 'LineStyle', '--', 'Color', [0.97, 0.55, 0.03]);
 end
 end
 

@Ternary_model/solve_tern_frap.m

@@ -3,8 +3,8 @@ function solve_tern_frap(T)
 % Assumes fixed profile of phi_tot
 fh_ic = @(x) flory_ic(x, T.a, T.u0, T.e, T.ic, T.x0);
 fh_bc = @(xl, ul, xr, ur, t) flory_bc(xl, ul, xr, ur, t, T.u0, T.e);
-fh_pde = @(x, t, u, dudx) flory_hugg_pde(x, t, u, dudx, T.a, T.b, T.e, T.u0,...
-                                         T.e_g0, T.u_g0);
+fh_pde = @(x, t, u, dudx) flory_hugg_pde(x, t, u, dudx, T.a, T.b, T.e,...
+                                          T.u0, T.e_g0, T.u_g0);
 T.sol = pdepe(T.geometry, fh_pde, fh_ic, fh_bc, T.x, T.t);
 end
 
@@ -38,7 +38,7 @@ elseif strcmp(ic, 'Gauss')
                 (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)));
     if x > a
-        u = p_out(D_p, D_m, ga, x0, x, 0.01/100);
+        u = p_out(D_p, D_m, ga, x0, x, 0.01/1000);
     else
         u = 0;
     end
@@ -51,9 +51,9 @@ function [pl,ql,pr,qr] = flory_bc(xl, ul, xr, ur, t, u0, e)
     pl = 0;
     ql = 1;
 %     % No flux
-%     pr = 0;%ur - 0.01;
-%     qr = 1;
+    pr = 0;%ur - 0.01;
+    qr = 1;
     % Dirichlet BC
-    pr = ur-u0+e;
-    qr = 0;
+%     pr = ur-u0+e;
+%     qr = 0;
 end
\ No newline at end of file

prob_laplace.m

@@ -85,38 +85,45 @@ end
 
 % different x0
 clear T;
-x0 = 0.01:0.01:10;
+x0 = 0.0:0.01:2.02;
 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} = 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}.solve_tern_frap()
 end
 toc
 %%
-ls = 0.01:0.01:2;
+ls = 0.0:0.0005:1;
 for i = 1:length(ls)
-q(i) = int_prob(ls(i), T, x0+6);
+q(i) = int_prob(ls(i), T, x0+5);
 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, p);
 plot(ls, q);
+%%
+int_prob(0.0, T, x0+5)
+%%
+for i = 1:7:length(T)
+    T{i}.plot_sim('plot', 1, 'red')
+end
 %% 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;
+    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 = 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
+            T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0)*p2;
+end
 end
-end
\ No newline at end of file