function output = runSubthresholdModel(obj, S, varargin)
%RUNSUBTHRESHOLDMODEL Run a simulation as defined by SubthresholdModel
%   O = runSubthresholdModel(PP, S) uses the output of subthrmodel to run 
%   new simulations mimicking that model as best as possible. The output O
%   contains six PricklyPear instances corresponding to the 6 DZW
%   conditions as well as all the spiketimes for each conditions.
%
%   runSubthresholdModel(..., 'axon_soma_distance', 45) sets the distance
%   between axon and soma (on ipsi dendrite) in microns. Default: 45

    % TODO Nicely implement this quickfix
    add_delay = 200; % ms

    % Handle input
    p = inputParser();
    p.addParamValue('axon_soma_distance', 45);
    p.parse(varargin{:});
    p = p.Results;

    % Handle output
    output = struct;
    output.irec = [];
    output.spt = struct;
    
    % The object must be an experiment (not a run)
    assert(~isempty(experiment(obj)), 'Please specify an experiment');
    assert(isempty(irec(obj)), 'Do not specify a run (irec)');
    
    % Log
    fprintf('Starting runSubthresholdModel...\n');
    
    % For each of the conditions
    for iCond = 1:numel(S.AllCond)
        % Get the cond
        cond = S.AllCond{iCond};
        
        % Log
        fprintf('Condition %i: %s\n', iCond, cond);
        fprintf('  - Setting parameters...\n');
        
        % Set PP parameters to S parameters
        obj.axon_soma_distance = p.axon_soma_distance;
        obj.eventtimes_contra = S.ev.(cond).C.evt + add_delay;
        obj.eventtimes_ipsi   = S.ev.(cond).I.evt + add_delay;
        obj.tstop = S.Param.Dur;
        obj.user_data = struct;
        obj.user_data.subthr = S;
        obj.user_data.cond = cond;
        
        % Log
        fprintf('  - Parameters set!\n');
        fprintf('  - Simulation duration: %i seconds.\n', ceil(obj.tstop/1e3));
        fprintf('  - Expected processing time: %i seconds.\n', ceil(obj.tstop/1e3)*5);
        fprintf('  - Running simulation...\n');
        
        % Run simulation
        run(obj);
        
        % Log
        fprintf('  - Run finished!\n');
        fprintf('  - Saving run...\n');
        
        % Save run
        nextIrec = save(child(obj, 0));
        output.irec(iCond) = nextIrec;
        
        % Log
        fprintf('  - Run saved as run #%i!\n', nextIrec);
        fprintf('  - Putting spiketimes in output...\n');
        
        % Put spiketimes in output struct
        output.spt.(cond) = struct;
        output.spt.(cond).spt = spiketimes(child(obj, nextIrec));
        
        % Log
        fprintf('  - Done!\n');
        
    end
    
    % Log
    fprintf('runSubthresholdModel is done!\n');

end