MATLAB code used to generate sixty different populations of particle arrangements and orientations

Description

clear all
clc
%**************************************************************************
% ****************Fracture toughness of different particles arrangements **************
% ****************************and orientations ********************************
% *************************************************************************
%**************************************************************************
%**************************************************************************
% Scriptfile name: Partarrang.m *************************************************
%**************************************************************************
% Puropse:
% This file is used to calculate the fracture toughness of different particles
% arrangement and orientations in front of crack tip of single edge notch
%
%**************************************************************************
PartArea=2513.275;
PartRadius=20;
ZoneVolume=9650972.63;
FirstIntegral_1=(PartArea*2)/ZoneVolume;
E=70000;
v=0.25;
%parameters
a = 0.00006;
b = 0.00002;
n=60;
x1_angle = 90; %REFERENCE ANGLE FOR THE X AXIS OF THE PARTICAL (relative to vertical);
%%%%%%%%%%%%%%%%%%%%%%%%%% Particle arrangements %%%%%%%%%%%%%%%%%%%%%%%%%%

x2_angle = 0; % 0 deg -> vertical particle
% 90 deg -> horizontal particle
% -45 deg -> diagonal right (positive slope)
% +45 deg -> diagonal left (negative slope)
for i=1:n
%%%%%%%%%%%%%%%%%%%%%%%%%% Particle (1)
PosArrang11(i) = randi(161)-1;%% O to 160
PosArrang21(i) = randi(81)-1;
%%%%%%%%%%%%%%%%%%%%%%%%%% Particle (2)
PosArrang12(i) = randi(161)-1;
PosArrang22(i) = randi(81)-1;
D = (((PosArrang12(i)-PosArrang11(i))^2 )+((PosArrang22(i)-
PosArrang21(i))^2))^0.5; %check
while D < 2*PartRadius
PosArrang11(i) = randi(161)-1;
PosArrang21(i) = randi(81)-1;
%%%%%%%%%%%%%%%%%%%%%%%%%% Particle (2)
PosArrang12(i) = randi(161)-1;
PosArrang22(i) = randi(81)-1;
D = (( (PosArrang12(i)-PosArrang11(i) )^2)+((PosArrang22(i)-
PosArrang21(i))^2))^0.5; %check
% Dist(i) = D;
end
while ( (PosArrang11(i) == 80 && PosArrang21(i) == 0) || (PosArrang12(i) ==
80 && PosArrang22(i) == 0))
PosArrang11(i) = randi(161)-1;%% O to 160
PosArrang21(i) = randi(81)-1;
%%%%%%%%%%%%%%%%%%%%%%%%%% Particle (2)
PosArrang12(i) = randi(161)-1;
PosArrang22(i) = randi(81)-1;
end
rArrang1(i)=sqrt((PosArrang11(i)-80)^2 + (PosArrang21(i))^2);
theta1(i)=atan2(PosArrang21(i),PosArrang11(i)-80);
rArrang2(i)=sqrt((PosArrang12(i)-80)^2 + (PosArrang22(i))^2);
theta2(i)=atan2(PosArrang22(i),PosArrang12(i)-80);
SecondIntegral1 = 1/((rArrang1(i))^(3/2))*cos(3*theta1(i)/2);
SecondIntegral2 = 1/((rArrang2(i))^(3/2))*cos(3*theta2(i)/2);
S1=abs( SecondIntegral1 * 1.28e-8);
S2=abs( SecondIntegral2 * 1.28e-8);

alpha1 = abs(x1_angle - theta1(i));
alpha2 = abs(x2_angle - theta2(i));
eps1 = 180/(70000*pi*a*b) * (pi*a *( a - (a+b)*cos(2*alpha1)));
eps2 = 180/(70000*pi*a*b) * (pi*a *( a - (a+b)*cos(2*alpha1)));
deltaK1(i)=((E*eps1)/(6*(sqrt(2*pi))*(1-
v)))*(FirstIntegral_1*0.000001)*S1;
deltaK2(i)=((E*eps2)/(6*(sqrt(2*pi))*(1-
v)))*(FirstIntegral_1*0.000001)*S2;
end