function [coords, L] = initCubicGrid(nPart,density)

        % Get the cooresponding box size
        L = (nPart/density)^(1/3);

        % Find the lowest perfect cube whose volume greater than or
        % equal to the number of particles

        nCube = ceil(nPart^(1/3));

        % Assign particle positions (Changed by Y.Zhao)
        xCoords = [0.5:nCube-0.5]*L/nCube;
        yCoords = xCoords;
        zCoords = xCoords;
        [X,Y,Z] = meshgrid(xCoords,yCoords,zCoords);
        X = X(:); X = X(1:nPart);
        Y = Y(:); Y = Y(1:nPart);
        Z = Z(:); Z = Z(1:nPart);
        coords = [X';Y';Z'];

    end

Published with MATLAB® R2012b