Random aggregate packing for concrete mesoscale analysis MATLAB

Given a 3D rectangular or cylindrical space with a certain volume fraction of aggregates (%) on the inside, I need to get them randomly distributed on the main volume according to a grading curve (which indicates the percentage of aggregates within the respective diameter ranges). The aggregates must be at least a distance d1 from the border of the domain and a distance d2 from each other (those values may change during the process). The aggregates must be regular octahedral shapes and I need the file to be exportable into a dxf file for further analysis in a Finite Element software (GiD). I already have the code implemented for 2d cases and for 3D using the take and place algorithm (https://www.researchgate.net/publication/223029627_Mesoscale_models_for_concrete_Homogenisation_and_damage_behaviour). It works for low density, however for higher aggregate densities, my code criteria has been poorly implemented (since Im comparing the distance between the spheres in which the aggregates are inscribed) and it is iterating for too long without completing the task. I have attached an illustration of the model as well as the zip folder with my current code, the main program is named "Rand_agg_generation_3D.m", the comments are in Portuguese.