clear
dx =5e-6; % x方向分辨率
dy =5e-6; % y方向分辨率
PML_size=20;%设置完美匹配层
% 阵列参数
element_size = [12, 12]; % x/y方向各6网格点，z方向1点
element_spacing = [1, 1]; % x/y间距1点，z方向无间隔
array_dims = [80, 1]; % x/y各80阵元，单层
% 总网格数
Nx=array_dims(1)*(element_size(1)+element_spacing(1)) - element_spacing(1)+2*PML_size;%x方向总网格数量
Ny=Nx;%y方向总网格数量
kgrid=kWaveGrid(Nx, dx, Ny, dy);
%设置介质材料属性
sound_speed_glass = 5600; % 玻璃声速 [m/s]
density_glass = 2400; % 玻璃密度 [kg/m^3]
sound_speed_water = 1500; % 水声速
density_water = 1000; % 水密度
sound_speed_air =343; % 空气声速
density_air = 1; % 空气密度
medium.sound_speed = sound_speed_glass*ones(Nx,Ny);%声速初始化
medium.density = density_glass*ones(Nx,Ny); %密度初始化
reflector_y_pos = 1+(PML_size + round(1e-3 / kgrid.dy));%反射面设置为在21个，避免反射波对计算区域造成影响
% 定义交替层的位置
x_start = PML_size + 1; %定义开始位置
x_end = Nx-PML_size; %定义结束位置
interval_num = 20; %分成20个区域
interval_grid_points = floor((Nx-2*PML_size+1) / interval_num); % 计算每个区域的网格点数，并使用floor函数取整
remaining_points =Nx-2*PML_size - interval_num * interval_grid_points; % 剩余的网格点数
target_column =x_start;
interval_num=20;
i = zeros(interval_num,2); %网格单元的长度
for interval_idx = 0:(interval_num-1)
current_x_start = x_start + interval_idx * interval_grid_points;
current_x_end = current_x_start + interval_grid_points - 1;
i(interval_idx+1,1)=current_x_start;
i(interval_idx+1,2)=current_x_end;
% 最后一层分配剩余点数
if interval_idx == interval_num - 1
current_x_end = current_x_end + remaining_points;
end
% 边界保护
if current_x_end > x_end
current_x_end = x_end;
end
% 交替材料设置
if mod(interval_idx, 2) == 0
medium.sound_speed(current_x_start:current_x_end, :, reflector_y_pos:end) = sound_speed_water;
medium.density(current_x_start:current_x_end, :, reflector_y_pos:end) = density_water;
else
medium.sound_speed(current_x_start:current_x_end, :, reflector_y_pos:end) = sound_speed_air;
medium.density(current_x_start:current_x_end, :, reflector_y_pos:end) = density_air;
end
end
imagesc(medium.sound_speed),
%设置换能器的位置和检测传感器的位置
tr.element_width = 12; % 换能器宽度所占网格数
tr.element_length= 1;% 换能器长度所占网格数
tr.element_spacing =1; % 换能器间隔所占网格数
tr.number_elements= 80; % 换能器数量
tr.radius=Inf;
transducer_width = tr.number_elements*tr.element_width+(tr.number_elements - 1) * tr.element_spacing;
tr.position = round([Nx/2 - transducer_width/2,Ny-PML_size]);
tr.active_elements=zeros(tr.number_elements,1);
tr.active_elements(1:80)=1;
% 定义激励超声源信号
%transducer.mask = transducer_mask;
source_strength = 1e6; % [Pa]
tone_burst_freq = 10e6; % [Hz]
tone_burst_cycles = 5;
t_end = 1e-6; % 设置仿真结束时间为10个周期
kgrid.makeTime(medium.sound_speed, [], t_end); % 使用kgrid的makeTime方法生成时间数组
% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
%input_signal=(source_strength./(medium.sound_speed* medium.density)).*input_signal;
tr.input_signal=input_signal;
transducer=kWaveTransducer(kgrid, tr);
sensor_data= kspaceFirstOrder2D(kgrid, medium, transducer,transducer);