Introducing plot, cleaning up non object-oriented file.

@Ternary_model/Ternary_model.m

 classdef Ternary_model < handle
+    % TERNARY_MODEL Numerical solution of ternary FRAP model with solvent, 
+    % bleached and unbleached species. Model is assumed to be equilibrated
+    % (bleached+unbleached=const.=pt). Then bleached species initial
+    % conditions are introduced. Integration of model via pdepe.
+
     properties
         pa = '/Users/hubatsch/ownCloud/Dropbox_Lars_Christoph/DropletFRAP/FRAP_paper/';
         g0
@@ -25,8 +30,11 @@ classdef Ternary_model < handle
                 T.e_g0 = params(5);
                 T.o_g0 = params(6);
             end
+            T.create_mesh();
+            T.phi_t = Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0);
         end
         create_mesh(T)
+        plot_sim(T)
         solve_tern_frap(T)
         spacing(T)
     end

@Ternary_model/plot_sim.m

+function plot_sim(T)
+%PLOT_SIM Show movie of simulation
+%% Plotting
+figure(1); hold on;
+for i = 1:length(T.t)
+    cla; xlim([-T.a-1, -T.a+3]); ylim([0, 0.7]); 
+    ax = gca;
+    ax.FontSize = 16;
+    xlabel('position'); ylabel('volume fraction');
+    plot(T.x, T.phi_tot(T.x, T.a, T.b, T.e, T.u0), 'LineWidth', 2, 'LineStyle', '--'); 
+    plot(T.x, T.sol(i, :), 'LineWidth', 2);
+    shg
+    pause();
+%     print([num2str(i),'.png'],'-dpng')
+end
+end
+

ternary_frap.m

-% Numerical solution of ternary FRAP model with solvent, bleached and
-% unbleached species. Model is assumed to be equilibrated
-% (bleached+unbleached=const.=pt). Then bleached species initial
-% conditions are introduced. Integration of model via pdepe.
-pa = '/Users/hubatsch/ownCloud/Dropbox_Lars_Christoph/DropletFRAP/FRAP_paper/';
-a = -1;
-b = 0.000025;
-u0 = 0.05;
-e = 0.4;
-e_g0 = 0.16;
-o_g0 = 0.2;
-%%
-g0 = gamma0(x, a, b, e_g0, o_g0);
-pt = phi_tot(x, a, b, e, u0);
-%% Make useful mesh (by inverting the tanh profile and using this as spacing)
-x = linspace(-a-1, -a+1, 3000);
-g = gamma0(x, a, 1*b, e_g0, o_g0);
-g_unique = unique(g);
-x = linspace(g_unique(1), g_unique(end-1), 300);
-g_inv = spacing(x, a, 2*b, e_g0, o_g0);
-g_inv = g_inv(2:end-1);
-x = [linspace(-a-1, g_inv(2), 30), g_inv(3:end-2), ...
-    linspace(g_inv(end-1), -a+1, 30), linspace(-a+1.1, 300, 3000)];
-%% Solve pde
-tic
-t = linspace(0, 2, 200);
-fh_ic = @(x) flory_ic(x, a, u0);
-fh_bc = @(xl, ul, xr, ur, t) flory_bc(xl, ul, xr, ur, t, u0);
-fh_pde = @(x, t, u, dudx) flory_hugg_pde(x, t, u, dudx, a, b, e, u0,...
-                                         e_g0, o_g0);
-sol = pdepe(0, fh_pde, fh_ic, fh_bc, x, t);
-toc
-%% Plotting
-figure(1); hold on;
-for i = 1:length(t)
-    cla; xlim([-a-1, -a+3]); ylim([0, 0.7]); 
-    ax = gca;
-    ax.FontSize = 16;
-    xlabel('position'); ylabel('volume fraction');
-    plot(x, phi_tot(x, a, b, e, u0), 'LineWidth', 2, 'LineStyle', '--'); 
-    plot(x, sol(i, :), 'LineWidth', 2);
-    shg
-    pause();
-%     print([num2str(i),'.png'],'-dpng')
-end
 %% Boundary condition: Jump makes sense
-[~, ind] = min(abs(x+a)); % Find x position of boundary
-plot(max(sol(:, 1:ind-1), [], 2)./min(sol(:, ind+1:end), [], 2), 'LineWidth', 2)
+[~, ind] = min(abs(T.x+T.a)); % Find x position of boundary
+plot(max(T.sol(:, 1:ind-1), [], 2)./min(T.sol(:, ind+1:end), [], 2), 'LineWidth', 2)
 make_graph_pretty('time', 'ratio outside/inside', 'unbleached component');
 figure; % Shouldn't be symmetric, different SS for phi_u and phi_b??
-plot(min(u0+e-sol(:, 1:ind-1), [], 2)./max(u0-sol(:, ind+1:end), [], 2), 'LineWidth', 2)
+plot(min(T.u0+T.e-T.sol(:, 1:ind-1), [], 2)./max(T.u0-sol(:, ind+1:end), [], 2), 'LineWidth', 2)
 make_graph_pretty('time', 'ratio outside/inside', 'bleached component');
+%% Plot flux at boundary
+
 %% Figures
 % Time course
 figure; hold on;
@@ -86,67 +43,5 @@ plot(s(1:2:100, :)', 'Color', [0.12156   0.4666666   0.7058823], 'LineWidth', 1.
 % plot(v_fit(1:2:100, :)', '--', 'Color', [1.   0.49803   0.0549], 'LineWidth', 1.5);
 make_graph_pretty('$x [\mu m]$', 'intensity [a.u.]', '');
 print('Timecourse_model_only.png','-dpng')
-%% Function definitions for pde solver
-function [c, f ,s] = flory_hugg_pde(x, t, u, dudx, a, b, e, u0,...
-                                    e_g0, o_g0)
-% Solve with full ternary model. Analytical derivatives.
-% pt ... phi_tot
-% gra_a ... analytical gradient of phi_tot
-
-pt = phi_tot(x, a, b, e, u0);
-gra_a = gradient_analytical(x, a, b, e);
-g0 = gamma0(x, a, b, e_g0, o_g0);
-c = 1;
-f = g0*(1-pt)/pt*(pt*dudx-u*gra_a);
-s = 0;
-end
-
-function u = flory_ic(x, a, u0)
-    % FRAP initial condition
-    if x < -a; u = 0.0; else; u = u0; end
-    % Peak outside initial condition
-%     if (x < -a+0.2) && (x > -a+0.1); u = 10; else; u = 0; end
-%     u = 0.3;
-%     u = phi_tot(x, -50, 1);
-%%
-%% Frank's solution to the transfer/rate problem via Laplace transform
-x0 = 3;
-D_m = 0.11;%g0(1)*(1-u0-e)/(u0+e); % to make equal to ternary FRAP
-D_p = 0.342;%g0(end)*(1-u0)/u0;
-ga = 1/9;
-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)));
-if x > -a
-    u = p_out(D_p, D_m, ga, x0, x, 3/100);
-else
-    u = 0;
-end
-end
-
-function [pl,ql,pr,qr] = flory_bc(xl, ul, xr, ur, t, u0)
-    pl = 0;
-    ql = 1;
-%     % No flux
-%     pr = 0;%ur - 0.01;
-%     qr = 1;
-    % Dirichlet BC
-    pr = ur - u0;
-    qr = 0;
-end
-
-function g0 = gamma0(x, a, b, e_g0, o_g0)
-    g0 = e_g0*(tanh((x+a)/b)+1)/2+o_g0;
-end
-
-function sp = spacing(x, a, b, e_g0, o_g0)
-    sp = b*atanh(2/(e_g0)*(x-o_g0)-1)-a;
-end
 
-function p = phi_tot(x, a, b, e, u0)
-    p = e*(tanh(-(x+a)/b)+1)/2+u0;
-end
 
-function grad = gradient_analytical(x, a, b, e)
-    grad = -e*(1-tanh(-(x+a)/b).^2)/b*0.5;
-end