clc;
clear;
% parameters setting
cubic_x = 80;
cubic_y = 80;
cubic_z = 80;
fiber_number = 100;
fiber_section = 50;
fiber_radius = 2;
point = struct('x_coordinate', 0, 'y_coordinate', 0, 'z_coordinate', 0, 'main_angle', 0, 'xy_angle', 0, 'theta_new', 0, 'phi_new', 0, 'mu0_x', 0, 'mu0_y', 0, 'mu0_z', 0, 'mu_i_x', 0, 'mu_i_y', 0, 'mu_i_z', 0, 'x_rotate', 0, 'y_rotate', 0, 'z_rotate', 0);
fibers = cell(fiber_number, fiber_section);
mu_0 = cell(fiber_number, 1);
mu_i = cell(fiber_number, fiber_section);
mu_new = cell(fiber_number, fiber_section);
for i = 1 : 1: fiber_number
for j = 1 : 1: fiber_section
fibers{i, j} = point;
end
end
kappa_1 = 10;
kappa_2 = 100;
% generate random numbers yielded to p(牟)
theta_0_array = zeros(fiber_number, 1);
p_array = zeros(fiber_number, 1);
B = 1;
flag = 1;
i = 1;
while flag == 1
x = rand(1) * pi;
y = rand(1) * 0.4;
p = (B * sin(x)) / (4 * pi * ((1 + (B * B - 1) * cos(x) * cos(x) )) ^1.5);
p_array(i, 1)= p;
if y <= p
theta_0_array(i, 1)= x;
i = i + 1;
end
if i == fiber_number + 1
flag = 0;
end
end
% generate the initial coordinate which is uniformly distributed in the cubic window
for i = 1 : 1 : fiber_number
x0 = 1 + rand(1) * (cubic_x - 1); % lowbound + (upbound - lowbound) * rand(1)
y0 = 1 + rand(1) * (cubic_y - 1);
z0 = 1 + rand(1) * (cubic_z - 1);
fibers{i, 1}.x_coordinate = x0;
fibers{i, 1}.y_coordinate = y0 / 4;
fibers{i, 1}.z_coordinate = z0;
fibers{i, 1}.main_angle = theta_0_array(i, 1);
fibers{i, 1}.xy_angle = 2 * pi * rand(1);
fibers{i, 1}.theta_new = 0;
fibers{i, 1}.phi_new = 0;
fibers{i, 1}.mu0_x = sin(fibers{i, 1}.main_angle) * cos(fibers{i, 1}.xy_angle);
fibers{i, 1}.mu0_y = sin(fibers{i, 1}.main_angle) * sin(fibers{i, 1}.xy_angle);
fibers{i, 1}.mu0_z = cos(fibers{i, 1}.main_angle);
mu_0{i, 1} = [fibers{i, 1}.mu0_x; fibers{i, 1}.mu0_y; fibers{i, 1}.mu0_z];
fibers{i, 1}.mu_i_x = 0;
fibers{i, 1}.mu_i_y = 0;
fibers{i, 1}.mu_i_z = 0;
mu_i{i, 1} = [0; 0; 0];
end
% generate the first new angle of the first new point with Fisher distribution
for i = 1 : 1 : fiber_number
R1 = rand(1);
R2 = rand(1);
kappa = sqrt((kappa_1 * fibers{i, 1}.mu0_x + kappa_2 * fibers{i, 1}.mu0_x)^2 + (kappa_1 * fibers{i, 1}.mu0_y + kappa_2 * fibers{i, 1}.mu0_y)^2 + (kappa_1 * fibers{i, 1}.mu0_z + kappa_2 * fibers{i, 1}.mu0_z)^2);
lambda = exp((-2) * kappa);
main_angle_temp = 2 * asin(sqrt((-log(R1 * (1 - lambda) + lambda)) / (2 * kappa)));
xy_angle_temp = 2 * pi * R2;
fibers{i, 2}.theta_new = acos(sin(main_angle_temp) * cos(xy_angle_temp) * fibers{i, 1}.mu0_x + sin(main_angle_temp) * sin(xy_angle_temp) * fibers{i, 1}.mu0_y + cos(main_angle_temp) * fibers{i, 1}.mu0_z);
fibers{i, 2}.phi_new = atan((cos(fibers{i, 1}.xy_angle) * sin(main_angle_temp) * sin(xy_angle_temp) - sin(fibers{i, 1}.xy_angle) * sin(main_angle_temp) * cos(xy_angle_temp)) / (cos(fibers{i, 1}.main_angle) * cos(fibers{i, 1}.xy_angle) * sin(main_angle_temp) * cos(xy_angle_temp) + cos(fibers{i, 1}.main_angle) * sin(fibers{i, 1}.xy_angle) * sin(main_angle_temp) * sin(xy_angle_temp) - sin(fibers{i, 1}.main_angle) * sin(main_angle_temp)));
fibers{i, 2}.main_angle = fibers{i, 1}.main_angle;
fibers{i, 2}.xy_angle = fibers{i, 1}.xy_angle;
fibers{i, 2}.mu0_x = fibers{i, 1}.mu0_x;
fibers{i, 2}.mu0_y = fibers{i, 1}.mu0_y;
fibers{i, 2}.mu0_z = fibers{i, 1}.mu0_z;
mu_0{i, 2} = [fibers{i, 2}.mu0_x; fibers{i, 2}.mu0_y; fibers{i, 2}.mu0_z];
fibers{i, 2}.mu_i_x = fibers{i, 2}.mu0_x;
fibers{i, 2}.mu_i_y = fibers{i, 2}.mu0_y;
fibers{i, 2}.mu_i_z = fibers{i, 2}.mu0_z;
mu_i{i, 2} = [fibers{i, 2}.mu_i_x; fibers{i, 2}.mu_i_y; fibers{i, 2}.mu_i_z];
fibers{i, 2}.x_coordinate = fibers{i, 1}.x_coordinate + 0.5 * fiber_radius * sin(fibers{i, 2}.theta_new) * cos(fibers{i, 2}.phi_new);
fibers{i, 2}.y_coordinate = fibers{i, 1}.y_coordinate + 0.5 * fiber_radius * sin(fibers{i, 2}.theta_new) * sin(fibers{i, 2}.phi_new);
fibers{i, 2}.z_coordinate = fibers{i, 1}.z_coordinate + 0.5 * fiber_radius * cos(fibers{i, 2}.theta_new);
end
% generate subsequent points from the second point to the last
for num = 1 : 1 : fiber_number
for i = 2 : 1 : (fiber_section - 1)
fibers{num, (i+1)}.main_angle = fibers{num, 1}.main_angle;
fibers{num, (i+1)}.xy_angle = fibers{num, 1}.xy_angle;
fibers{num, (i+1)}.mu0_x = fibers{num, 1}.mu0_x;
fibers{num, (i+1)}.mu0_y = fibers{num, 1}.mu0_y;
fibers{num, (i+1)}.mu0_z = fibers{num, 1}.mu0_z;
mu_0{num, (i+1)} = [fibers{num, (i+1)}.mu0_x; fibers{num, (i+1)}.mu0_y; fibers{num, (i+1)}.mu0_z];
x_temp_main = sin(fibers{num, 1}.main_angle) * cos(fibers{num, 1}.xy_angle);
y_temp_main = sin(fibers{num, 1}.main_angle) * sin(fibers{num, 1}.xy_angle);
z_temp_main = cos(fibers{num, 1}.main_angle);
x_temp_last = sin(fibers{num, (i-1)}.theta_new) * cos(fibers{num, (i-1)}.phi_new);
y_temp_last = sin(fibers{num, (i-1)}.theta_new) * sin(fibers{num, (i-1)}.phi_new);
z_temp_last = cos(fibers{num, (i-1)}.theta_new);
R1 = rand(1);
R2 = rand(1);
kappa_new = sqrt((kappa_1 * x_temp_main + kappa_2 * x_temp_last)^2 + (kappa_1 * y_temp_main + kappa_2 * y_temp_last)^2 + (kappa_1 * z_temp_main + kappa_2 * z_temp_last)^2);
lambda_new = exp((-2) * kappa_new);
main_angle_temp_new = 2 * asin(sqrt((-log(R1 * (1 - lambda_new) + lambda_new)) / (2 * kappa_new)));
xy_angle_temp_new = 2 * pi * R2;
fibers{num, (i+1)}.theta_new = acos(sin(main_angle_temp_new) * cos(xy_angle_temp_new) * x_temp_main + sin(main_angle_temp_new) * sin(xy_angle_temp_new) * y_temp_main + cos(main_angle_temp_new) * z_temp_main);
fibers{num, (i+1)}.phi_new = atan((cos(fibers{num, 1}.xy_angle) * sin(main_angle_temp_new) * sin(xy_angle_temp_new) - sin(fibers{num, 1}.xy_angle) * sin(main_angle_temp_new) * cos(xy_angle_temp_new)) / (cos(fibers{num, 1}.main_angle) * cos(fibers{num, 1}.xy_angle) * sin(main_angle_temp_new) * cos(xy_angle_temp_new) + cos(fibers{num, 1}.main_angle) * sin(fibers{num, 1}.xy_angle) * sin(main_angle_temp_new) * sin(xy_angle_temp_new) - sin(fibers{num, 1}.main_angle) * cos(main_angle_temp_new)));
fibers{num, (i+1)}.mu_i_x = sin(fibers{num, (i+1)}.main_angle) * cos(fibers{num, (i+1)}.xy_angle);
fibers{num, (i+1)}.mu_i_y = sin(fibers{num, (i+1)}.main_angle) * sin(fibers{num, (i+1)}.xy_angle);
fibers{num, (i+1)}.mu_i_z = cos(fibers{num, (i+1)}.main_angle);
mu_i{num, (i+1)} = [fibers{num, (i+1)}.mu_i_x; fibers{num, (i+1)}.mu_i_y; fibers{num, (i+1)}.mu_i_z];
fibers{num, (i+1)}.x_coordinate = fibers{num, i}.x_coordinate + 0.5 * fiber_radius * sin(fibers{num, (i+1)}.theta_new) * cos(fibers{num, (i+1)}.phi_new);
fibers{num, (i+1)}.y_coordinate = fibers{num, i}.y_coordinate + 0.5 * fiber_radius * sin(fibers{num, (i+1)}.theta_new) * sin(fibers{num, (i+1)}.phi_new);
fibers{num, (i+1)}.z_coordinate = fibers{num, i}.z_coordinate + 0.5 * fiber_radius * cos(fibers{num, (i+1)}.theta_new);
if ((fibers{num, (i+1)}.x_coordinate > cubic_x) || (fibers{num, (i+1)}.x_coordinate < 0) || (fibers{num, (i+1)}.y_coordinate > cubic_y) || (fibers{num, (i+1)}.y_coordinate < 0) || (fibers{num, (i+1)}.z_coordinate > cubic_z) || (fibers{num, (i+1)}.z_coordinate < 0))
fibers{num, (i+1)}.x_coordinate = fibers{num, i}.x_coordinate;
fibers{num, (i+1)}.y_coordinate = fibers{num, i}.y_coordinate;
fibers{num, (i+1)}.z_coordinate = fibers{num, i}.z_coordinate;
continue
end
end
end
% calculate the mean fiber orientation of a created ball chain
n = cell(fiber_number, 1);
alpha_angle= cell(fiber_number, 1);
fo
