solve_tern_frap.m

function solve_tern_frap(T)
%SOLVE_TERN_FRAP solves Fokker Planck equation for ternary FRAP
% Assumes fixed profile of phi_tot
tic
fh_ic = @(x) flory_ic(x, T.a, T.b, T.u0, T.e, T.ic, T.x0, T.ic_c);
fh_bc = @(xl, ul, xr, ur, t) flory_bc(xl, ul, xr, ur, t, T.u0, T.e, T.j);
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.v, T.j,...
                                          T.mode);
opts = odeset('RelTol',1e-3,'AbsTol',1e-6);
T.sol = pdepe(T.geometry, fh_pde, fh_ic, fh_bc, T.x, T.t, opts);
toc
end

%% Function definitions for pde solver
function u = flory_ic(x, a, b, u0, e, ic, x0, ic_c)
if strcmp(ic, 'FRAP')
    % FRAP initial condition
    if x < -a
        u = ic_c*Ternary_model.phi_tot(x, a, b, e, u0, 0); % normalize
    else
        % Concentration outside, away from droplet
        u = u0-e;
%         u = Ternary_model.phi_tot(x, a, b, e, u0, 0);
    end
elseif strcmp(ic, 'Gauss')
    % 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))./...
                (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/1000);
    else
        u = 0;
    end
elseif strcmp(ic, 'phi_tot')
    % This uses standard phi_tot profile, which is not correct for moving
    % boundary with client
%     u = Ternary_model.phi_tot(x, a, 0.021282, 0.460977, u0, 0);

    % correct form of initial profile for moving boundary with client
    chi = -4.530864768482371;
    tern = @(x) Ternary_model.phi_tot(x, a, 0.021282, 0.460977, u0, 0);
    exp_tern = @(x) exp(-chi*Ternary_model.phi_tot(x, a, 0.021282, 0.460977, u0, 0));
    off = tern(200);
    P = tern(-200)/tern(200);
    A = (P-1)*off/(exp_tern(-200)-exp_tern(200));
    B = off - A*exp_tern(200);
    u =  A*exp_tern(x)+B;
elseif strcmp(ic, 'Substrate')
    u = zeros(2, 1);
    if x < -2*a
        u(1) = 0;
    else
        u(1) = 1;
    end
    u(2) = 0;
elseif isnumeric(ic)
    % use ic as initial condition by interpolating the given values between
    % 0 and a
    if x < -a
        u = interp1(linspace(0, -a, length(ic)), ic, x);
    else
        u = u0-e;
    end
else
    disp('Initial condition not supported');
end
end

function [pl,ql,pr,qr] = flory_bc(xl, ul, xr, ur, t, u0, e, j)
    if j == 0
        pl = 0;
        ql = 1;
        % No flux
        pr = 0;%ur - 0.01;
        qr = 1;
        % Dirichlet BC
%         pr = ur-u0+e;
%         qr = 0;
    elseif j == -1 %% For now this means we're looking at Hammerhead/Ribo
        pl = [0; 0];
        ql = [1; 1];
        % No flux
        pr = [0; 0];%ur - 0.01;
        qr = [1; 1];
    % Dirichlet BC constant flux case
    else
        pr = ur;
        qr = 0;
        pl = ul;
        ql = 0;
    end
end