Не могу написать программу на матлабе, для генерации кристаллической решетки сплава

Помогите написать код.Нужно создать кристаллическую решетку сплава металлов и случайное распределение атомов в металле.Кристаллическую решетку сгенерировал, а сделать перестановку в кристаллической решетки не получается. Код на Матлабе:

% Параметры моделирования
N=4000; %Количество атомов
L=992; %Размер решётки
steps = 100; %Количество шагов Монте-Карло
%Инициализация начальной конфигурации решётки
x = randi(10,50,50);
y = randi(10,50,50);
z = randi(10,50,50);
x(1,1,1)=1;
y(1,1,1)=1;
z(1,1,1)=1;
x(2,1,1)=10;
y(2,1,1)=1;
z(2,1,1)=1;
x(3,1,1)=1;
y(3,1,1)=10;
z(3,1,1)=1;
x(4,1,1)=10;
y(4,1,1)=10;
z(4,1,1)=1;
x(5,1,1)=1;
y(5,1,1)=1;
z(5,1,1)=10;
x(6,1,1)=10;
y(6,1,1)=1;
z(6,1,1)=10;
x(7,1,1)=1;
y(7,1,1)=10;
z(7,1,1)=10;
x(8,1,1)=10;
y(8,1,1)=10;
z(8,1,1)=10;
x(9,1,1)=10;
y(9,1,1)=1;
z(9,1,1)=1;
k = randi(8,9,9);
l = randi(8,9,9);
m = randi(8,9,9);
q=[];
w=[];
e=[];
for i = 1:9
for g = 1:9
q = [q, 0.5 + i];
w = [w, 0.5 + g];
e = [e, 1.5];
end
end

for i = 1:9
for g = 1:9
q = [q, 0.5 + i];
w = [w, 0.5 + g];
e = [e, 2.5];
end
end
for i = 1:9
for g = 1:9
q = [q, 0.5 + i];
w = [w, 0.5 + g];
e = [e, 3.5];
end
end
for i = 1:9
for g = 1:9
q = [q, 0.5 + i];
w = [w, 0.5 + g];
e = [e, 4.5];
end
end
for i = 1:9
for g = 1:9
q = [q, 0.5 + i];
w = [w, 0.5 + g];
e = [e, 5.5];
end
end
for i = 1:9
for g = 1:9
q = [q, 0.5 + i];
w = [w, 0.5 + g];
e = [e, 6.5];
end
end
for i = 1:9
for g = 1:9
q = [q, 0.5 + i];
w = [w, 0.5 + g];
e = [e, 7.5];
end
end
for i = 1:9
for g = 1:9
q = [q, 0.5 + i];
w = [w, 0.5 + g];
e = [e, 8.5];
end
end
for i = 1:9
for g = 1:9
q = [q, 0.5 + i];
w = [w, 0.5 + g];
e = [e, 9.5];
end
end
hold on
scatter3(x,y,z);
scatter3(q,w,e);
hold off
xlabel('x');
ylabel('y');
zlabel('z');
title ('Структура')

Посмотреть все изображения

Пробовал через методы монте-карло. но ничего хорошего не получилось. Мб у вас получится:

% Параметры моделирования
N = 200; % Количество атомов
L = 200; % Размер решётки
steps = 1; % Количество шагов Монте-Карло

% Инициализация начальной конфигурации решётки
x = zeros(L, L, L);
y = zeros(L, L, L);
z = zeros(L, L, L);
x(1, 1, 1) = 0;
y(1, 1, 1) = 0;
z(1, 1, 1) = 0;
x(2, 1, 1) = L;
y(2, 1, 1) = 0;
z(2, 1, 1) = 0;
x(3, 1, 1) = 0;
y(3, 1, 1) = L;
z(3, 1, 1) = 0;
x(4, 1, 1) = L;
y(4, 1, 1) = L;
z(4, 1, 1) = 0;
x(5, 1, 1) = 0;
y(5, 1, 1) = 0;
z(5, 1, 1) = L;
x(6, 1, 1) = L;
y(6, 1, 1) = 0;
z(6, 1, 1) = L;
x(7, 1, 1) = 0;
y(7, 1, 1) = L;
z(7, 1, 1) = L;
x(8, 1, 1) = L;
y(8, 1, 1) = L;
z(8, 1, 1) = L;
x(9, 1, 1) = 1;
y(9, 1, 1) = 0;
z(9, 1, 1) = 0;

for h = 1:L
    for q = 1:L
        for i = 10:N
            x(i, q, h) = x(9, 1, 1) + 1;
            y(i, q, h) = 0;
            z(i, q, h) = 0;
        end
    end
end

for i = 10:(L-1)
    for q = 1:(L-1)
        for h = 1:(L-1)
            if (q == 1) && (h == 1)
                x(i+1, q, h) = x(i, q, h) + 1;
                y(i+1, q, h) = 0;
                z(i+1, q, h) = 0;
            elseif (q > 1) && (h == 1)
                x(i+1, q, h) = x(i, q, h) + 1;
                y(i+1, q, h) = y(i, q, h) + 1;
                z(i+1, q, h) = 0;
            elseif (q > 1) && (h > 1)
                x(i+1, q, h) = x(i, q, h) + 1;
                y(i+1, q, h) = y(i, q, h) + 1;
                z(i+1, q, h) = z(i, q, h) + 1;
            end
        end
    end
end

for h = (L):L
    for q = (L):L
        for i = (L):L
            x(i, q, h) = 0.5;
            y(i, q, h) = 0.5;
            z(i, q, h) = 0;
        end
    end
end

for h = (L-1):(2*L-1)
    for q = (L):(2*L-1)
        for i = (L-1):(2*L-1)
            if (q == (L-1)) && (h == (L))
                x(i+1, q, h) = x(i, q, h) + 1;
                y(i+1, q, h) = 0;
                z(i+1, q, h) = 0;
            elseif (q > (L-1)) && (h == (L))
                x(i+1, q, h) = x(i, q, h) + 1;
                y(i+1, q, h) = y(i, q, h) + 1;
                z(i+1, q, h) = 0;
            elseif (q > (L-1)) && (h > (L))
                x(i+1, q, h) = x(i, q, h) + 1;
                y(i+1, q, h) = y(i, q, h) + 1;
                z(i+1, q, h) = z(i, q, h) + 1;
            end
        end
    end
end

% Разворачивание трехмерного массива в векторы
x = reshape(x, N, 1);
y = reshape(y, N, 1);
z = reshape(z, N, 1);

% Визуализация решётки
scatter3(x, y, z);
xlabel('x');
ylabel('y');
zlabel('z');
title('Структура');

Ошибка Матлаба:

Index in position 2 exceeds array bounds. Index must not exceed 200.
Error in untitled (line 86)
                x(i+1, q, h) = x(i, q, h) + 1;