code.rar_随机游走_随机游走模型

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