Hi All,
I would like to define 4 layers medium and a ball within the layers.
I could assign the properties of the layers, but how to assign the properties of the ball?
Here is my code:
% define a random distribution of scatterers for the medium
background_map_mean = 1;
background_map_std = 0.008;
% define the thicknesses of layer phantom in x direction
d_layer_1 = 5.4e-3; % [m]
d_layer_2 = 4.9e-3; % [m]
d_layer_3 = 4.4e-3; % [m]
d_layer_4 = 6.3e-3; % [m]
% define the thicknesses of layer phantom in x direction
position_x1 = 1:round(Nx*(d_layer_1/x));
position_x2 = round(Nx*(d_layer_1/x))+1:round(Nx*((d_layer_1+d_layer_2)/x));
position_x3 = round(Nx*((d_layer_1+d_layer_2)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x));
position_x4 = round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x));
% define the width of phantom in y and z direction
width_y = 29e-3; % width_phantom_in y direction [grid points]*dy = [m]
width_z = 35e-3; % width_phantom_in z direction [grid points]*dz = [m]
% define position of each layer in y and z direction
position_y = round(Ny*(((y-width_y)/2)/y))+1:round(Ny*(y-(y-width_y)/2)/y);
position_z = round(Nz*(((z-width_z)/2)/z))+1:round(Nz*(z-(z-width_z)/2)/z);
% define the dimension of each layer (equal to its background map)
width_backgroundmap_layer_1 = [round(Nx*(d_layer_1/x)), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))];
width_backgroundmap_layer_2 = [round(Nx*((d_layer_1+d_layer_2)/x))-(round(Nx*(d_layer_1/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))];
width_backgroundmap_layer_3 = [round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))-(round(Nx*((d_layer_1+d_layer_2)/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))];
width_backgroundmap_layer_4 = [round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x))-(round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))];
% % define the properties of the propagation medium
% % water
medium.sound_speed = c0*ones(Nx, Ny, Nz); % [m/s]
medium.density = rho0*ones(Nx, Ny, Nz); % [kg/m^3]
% layer_1
background_map_layer_1 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_1);
medium.sound_speed(position_x1, position_y, position_z) = 1520.*background_map_layer_1; % [m/s]
medium.density(position_x1, position_y, position_z) = 1281.*background_map_layer_1; % [kg/m^3]
% layer_2
background_map_layer_2 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_2);
medium.sound_speed(position_x2, position_y, position_z) = 1218.*background_map_layer_2; % [m/s]
medium.density(position_x2, position_y, position_z) = 1510.*background_map_layer_2; % [kg/m^3]
% layer_3
background_map_layer_3 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_3);
medium.sound_speed(position_x3, position_y, position_z) = 1325.*background_map_layer_3; % [m/s]
medium.density(position_x3, position_y, position_z) = 1643.*background_map_layer_3; % [kg/m^3]
% layer_4
background_map_layer_4 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_4);
medium.sound_speed(position_x4, position_y, position_z) = 1429.*background_map_layer_4; % [m/s]
medium.density(position_x4, position_y, position_z) = 1684.*background_map_layer_4; % [kg/m^3]
% define a ball for a highly scattering region
radius = 4e-3; % [m]
ball = makeBall(Nx, Ny, Nz, Nx/2, Ny/2, Nz/2, round(radius/dx));
% define properties of the ball
c_ball = 3000;
rho_ball = 2500;
% assign region of the ball
medium.sound_speed(ball) = c_ball; % [m/s]
medium.density(ball) = rho_ball; % [kg/m^3]I got an error message 'Subscript indices must either be real positive integers or logicals' for the last 2 lines.
Do you have a suggestion how to do it properly?
Many thanks.
Maktjik