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);
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, T.v,...
                                          T.mode);
T.sol = pdepe(T.geometry, fh_pde, fh_ic, fh_bc, T.x, T.t);
toc
end

%% Function definitions for pde solver


function u = flory_ic(x, a, b, u0, e, ic, x0)
if strcmp(ic, 'FRAP')
    % FRAP initial condition
    if x < -a
        u = 0.0;
    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')
    u = Ternary_model.phi_tot(x, a, 0.021282, 0.460977, u0, 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)
    pl = 0;
    ql = 1;
%     % No flux
    pr = 0;%ur - 0.01;
    qr = 1;
    % Dirichlet BC
%     pr = ur-u0+e;
%     qr = 0;
end