% b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
% e = @(chi) sqrt(3/8*(chi-2));
% t = linspace(0, 10, 100);
% T = Ternary_model(0, 'FRAP', {-0.1, b(7.9/3, 10^-12), 0.5, e(7.9/3),...
%                                      0, 1, 30, 7, 0, 'Constituent', 0},...
%                                      t, 2);
% T.solve_tern_frap()
load ~/Desktop/f_eos_parts_local.mat
fit_seed_part = fit_seed;
load ~/Desktop/f_eos.mat
%% Select individual simulation
jj = find(n_f==25 & fix_params(3, :)'==2);
jj = jj(20);
D = abs(f_temp{n_f(jj)}.conv_factor_D*f_temp{n_f(jj)}.x_seed);
[cost, T] = to_min(fit_seed(jj), fit, [D; fix_params(2:3, jj)], fixed,...
                   f_temp{n_f(jj)});
%% Analytical solution for steady state (see note from 03/11/20 in FRAP
D_out = (1-T.phi_t(end))*T.ga0(end);
L = sqrt(D_out)*80;
R = -T.a+10*T.b;
t_ind = 11;
ind_R = find(T.x>R, 1, 'first');
x = T.x(ind_R:end);
j = -D_out*(T.sol(t_ind, ind_R+1)-T.sol(t_ind, ind_R))/...
                      (T.x(ind_R+1)-T.x(ind_R));
B = j*R^2/D_out;
% infinite boundary
A = T.phi_t(end);
% finite boundary
% ind_L = find(T.x>L-T.a, 1);
% A = T.sol(t_ind, ind_L)-B/T.x(ind_L);
c = A+B./x;
T.plot_sim('plot', 10, 'green', 1);
plot(x, c, 'r', 'LineWidth', 3);
diff = sum(T.sol(t_ind, ind_R:end)-c)/length(c);
axis([-inf, inf, -inf, inf]);
j_qs = -D_out*(c(2)-c(1))/(x(2)-x(1));