function t_halves = FRAP_biexpfit(CellOfFRAPData, TimePoints, filename_for_writing)
    % Computes t_halves from normalized FRAP curves
    % Output - 
    %    t_halves = 1 x k array where k = number of FRAP curves in
    %    CellOfFRAPData
    %
    % Inputs - 
    %    CellOfFRAPData - cell of FRAP curves. Each element is a 1 x k array 
    %                     where k denotes the total number of frames 
    %    TimePoints - 0 to the total time of the simulation (in seconds)
    %    filename_for_writing - any suitable filename (with extension)
    % 
    % Example - 
    %    t_halves = FRAP_biexpfit({FRAPData1, FRAPData2], (0:6000)*0.00005,'t_halves_log.txt')
    
    ft_eqn = fittype('I0 - alpha*exp(-beta*x) - gamma*exp(-delta*x)'); 
    % Biexp fit equation
    syms f(x)
    t_halves = zeros(1,length(CellOfFRAPData));
    fileID = fopen(filename_for_writing,'a');
    for i = 1:length(CellOfFRAPData)
        [ft,gof] = fit(TimePoints',CellOfFRAPData{i}(4000:end)',ft_eqn,...
            'Algorithm','Levenberg-Marquardt');
        f(x) = ft.I0 - ft.alpha*exp(-ft.beta*x) - ft.gamma*exp(-ft.delta*x);
        try
            t_halves(i) = double(solve(f(x)==ft.I0/2,x));
            fprintf(fileID, "Entry %.1g\n",i);
            fprintf(fileID, 'R-square value is: %.3g\n', gof.rsquare);
            fprintf(fileID, "t1/2 = %.3g s\n************\n", ...
                t_halves(i)); % 
        catch
            continue
        end % try catch is needed as sometimes the fit is wrong
        % Call this function multiple times, evaluate
    end
    fclose(fileID);
end