function [cost, v_fit] = to_min_theory(r_n, t, x_fit, reference)
% TO_MIN_THEORY Takes reference simulation and tries to find simulation 
% that fits total intensity recovery. Uses the full model, including bulk.
    T = Ternary_model(2, 'FRAP', [-1, 0.001, 0.5, 0.4, x_fit(1), x_fit(2), 40], t, 1);
    T.solve_tern_frap();
    v_fit = T.sol(:, 1:100);
% Return cost as the sum of squared differences of model and experiment
% cost = sum(sum((v_fit-r_n).^2));
cost = sum((sum(v_fit, 2) - reference).^2);
end