function wh = submit_simulations(wh,options)

arguments
    wh 
    options.parallel = 1;
end

%%% 2) SUBMIT SIMULATION
% Initialize job results
if options.parallel
    if isempty(gcp('nocreate'))
        parpool;
    end
    results = parallel.FevalFuture.empty();
else
    results = [];
end


lin_idx = 1;
for lin_idx = 1:wh.getLastLinIndice
    optionalVars = struct();
    if ~isempty(wh.getDimension)
        % Build the optionalVars struct
        [parametervalues,parameternames]=wh.getPhysIndicesByLinIndex(lin_idx);
        for pidx = 1:numel(parameternames)
            optionalVars.(parameternames{pidx}) = parametervalues{pidx};
        end
    end

    %%% SIMULATION HERE
    if options.parallel

        numOutputs = 1;
        results(lin_idx) = parfeval(@imdd_simulation_minimal, numOutputs, optionalVars);

    else

        finalresults{lin_idx} = imdd_simulation_minimal(optionalVars);
        wh.addValueToStorageByLinIdx(finalresults{lin_idx}, 'ber', lin_idx);

    end
end



if options.parallel

    %%% 4) Setup waitbar
    h = waitbar(0, 'Processing Simulations...');
    % Helper function to compute progress
    updateWaitbar = @(~) waitbar(mean(arrayfun(@(f) strcmp(f.State, 'finished'), results)), h);

    fprintf('Fetching results... \n');

    % Update the waitbar after each simulation
    updateWaitbarFutures = afterEach(results, updateWaitbar, 0);

    % Close the waitbar after all simulations complete
    afterAll(updateWaitbarFutures, @(~) delete(h), 0);

    %%% 7) Fetch final results after all computations
    fetchOutputs(results);

    for ridx = 1:length(results)
        wh.addValueToStorageByLinIdx(results(ridx).OutputArguments{1}, 'ber', ridx);
    end

end


end