Meshing now symmetric around boundary. Gamma0=1 now implicit.

@Ternary_model/create_mesh.m

@@ -6,27 +6,27 @@ function create_mesh(T)
 % x = [0, 0.01/T.precision/10];%/10
 x = linspace(0, 0.75, 40);
 x = [x, x(end)+0.01/T.precision/10];%/10
-while x(end)<-2*T.a
+while x(end)<-T.a
     x_temp = x(end)-x(end-1);
     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.002/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.001/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
+    while phi_nm1-phi_n > 0.002/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
     x = [x, x(end)+x_temp];
 end
+x = [x, x(end)+cumsum(fliplr(diff(x)))];
 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);
 T.x = [x, x_middle, x_right];
 T.x = T.x(1:find(T.x-T.system_size>=0, 1)-1);
+
+% Fine mesh close to x0 if ic = 'Gauss'
+[~, ind] = min(abs(T.x-T.x0));
+if T.x(ind+20)-T.x(ind) > 0.1
+    T.x = [T.x(1:ind-21), linspace(T.x(ind-20), T.x(ind+20),...
+                                        ceil((T.x(ind+20)-T.x(ind-20))*200)),...
+                          T.x(ind+21:end)];
 end
 

@Ternary_model/plot_sim.m

@@ -13,7 +13,7 @@ if strcmp(mode, 'mov')
         cla; xlim([-T.a-1, -T.a+3]); ylim([0, 0.2]);
         xlim([0, 100]);
         ax = gca;
-        ax.FontSize = 16;
+        ax.FontSize = 12;
         xlabel('position'); ylabel('probability [not normalized]');
         plot(T.x, T.phi_tot(T.x, T.a, T.b, T.e, T.u0), 'LineWidth', 2, ...
             'LineStyle', '--');
@@ -26,9 +26,9 @@ if strcmp(mode, 'mov')
 elseif strcmp(mode, 'plot')
     figure(20);
     hold on;
-    xlim([-inf, 3]); ylim([0, inf]);
+    xlim([-inf, 10]); ylim([-inf, 2]);
     ax = gca;
-    ax.FontSize = 26;
+    ax.FontSize = 12;
     xlabel('x [\mum]'); ylabel('volume fraction'); 
     plot(T.x, T.sol([1, 2:nth:si(1)], :), 'LineWidth', 1, 'Color', color);
 %     plot(T.x, T.sol(nth, :), 'LineWidth', 1, 'Color', color);

@Ternary_model/solve_tern_frap.m

@@ -18,7 +18,6 @@ function [c, f ,s] = flory_hugg_pde(x, t, u, dudx, a, b, e, u0,...
 pt = Ternary_model.phi_tot(x, a, b, e, u0);
 gra_a = Ternary_model.gradient_analytical(x, a, b, e);
 g0 = Ternary_model.gamma0(x, a, b, e_g0, u_g0);
-% g0 = 1;
 c = 1;
 f = g0*(1-pt)/pt*(pt*dudx-u*gra_a);
 s = 0;
@@ -34,10 +33,9 @@ if strcmp(ic, 'FRAP')
         u = u0-e;
     end
 elseif strcmp(ic, 'Gauss')
-    % Frank/Jonathan/Stefano solution to the transfer/rate problem via 
-    % Laplace transform
-    D_m = 1*(1-u0-e); % to make equal to ternary FRAP
-    D_p = 1*(1-u0+e);
+    % Frank/Jonathan/Stefano solution via Laplace transform
+    D_m = (1-u0-e); % to make equal to ternary FRAP
+    D_p = (1-u0+e);
     ga = (u0-e)/(u0+e);
     p_out = @(D_p, D_m, ga, x0, x, t) 1./(sqrt(4*D_p*pi*t))*...
                 (exp(-(x+x0).^2./(4*D_p*t))*(ga*sqrt(D_p)-sqrt(D_m))./...