Hi, I've defined a 3d grid and trying to run the kspaceFirstOrder3D function. However, the error message shows "Error using expandMatrix
exp_coeff must be a 1 or 2 element array.
Error in kspaceFirstOrder_expandGridMatrices (line 160)
source.p0 = expandMatrix(source.p0, expand_size, 0);
Error in kspaceFirstOrder_inputChecking (line 1558)
kspaceFirstOrder_expandGridMatrices;
Error in kspaceFirstOrder3D (line 574)
kspaceFirstOrder_inputChecking;
Error in kWave_PSI_simulation (line 102)
sensorDataTmp = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});"
I've checked the example script as well, but same issue is happening. It seems that this is something related to built in functions.. I appreciate any solutions. Here is my script for your reference.
clc; clear; close all;
%%
c = 1540;
fs = 30e6;
% fc = 2e6;
% band_width = 50; % percent
% pass_band = nFc * (1+[-band_width band_width]/100) / nFs * 2;
bRecord = false;
% computation settings
file_name = 'PSI_simulation';
input_args = {'DataCast', 'single',...
'PlotLayout',true,...
'PMLInside',false,...
'PlotPML',false,...
'RecordMovie', bRecord,...
'MovieName', file_name,...
'MovieProfile', 'MPEG-4'};
%%
dir_cur = pwd;
dirSave = [dir_cur '/../data/Acoustic_simulation/' file_name];
mkdir(dirSave);
%% dimension & phantom setup
% size of the computational grid
grid_size = 2*[40 40 10]*1e-3; % x, y, z
grid_interval = [1 1 1]*5e-4;
Nx = round(grid_size(1)/grid_interval(1));
Ny = round(grid_size(2)/grid_interval(2));
Nz = round(grid_size(3)/grid_interval(3));
kgrid = kWaveGrid(Nx, grid_interval(1), Ny, grid_interval(2), Nz, grid_interval(3));
axis_x = linspace(-0.5*(Nx-1), 0.5*(Nx-1), Nx)*grid_interval(1);
axis_y = linspace(-0.5*(Ny-1), 0.5*(Ny-1), Ny)*grid_interval(2);
axis_z = linspace(0, Nz-1, Nz)*grid_interval(3);
%% scanning sequence
scanning_interval = [1 1]*1e-3; % x, y
scanning_pos_x = -0.5*(grid_size(1)/2):scanning_interval(1):0.5*(grid_size(1)/2);
scanning_pos_y = -0.5*(grid_size(2)/2):scanning_interval(2):0.5*(grid_size(2)/2);
[scanning_pos_y, scanning_pos_x] = ndgrid(scanning_pos_x, scanning_pos_y);
scanning_pos_y_grid = round((scanning_pos_y + grid_size(2)/2)/grid_interval(2))+1;
scanning_pos_x_grid = round((scanning_pos_x + grid_size(1)/2)/grid_interval(1))+1;
%% target
target_pos = [0 0 1; 0 10 1; 10 0 1; 10 10 1; 0 -10 1; -10 0 1; -10 -10 1; -10 10 1; 10 -10 1;
% 0 0 2; 0 1 2; 1 0 2; 1 1 2; 0 -1 2; -1 0 2; -1 -1 2; -1 1 2; 1 -1 2;
0 0 3; 0 10 3; 10 0 3; 10 10 3; 0 -10 3; -10 0 3; -10 -10 3; -10 10 3; 10 -10 3;
% 0 0 4; 0 1 4; 1 0 4; 1 1 4; 0 -1 4; -1 0 4; -1 -1 4; -1 1 4; 1 -1 4;
0 0 5; 0 10 5; 10 0 5; 10 10 5; 0 -10 5; -10 0 5; -10 -10 5; -10 10 5; 10 -10 5;
% 0 0 6; 0 1 6; 1 0 6; 1 1 6; 0 -1 6; -1 0 6; -1 -1 6; -1 1 6; 1 -1 6;
0 0 7; 0 10 7; 10 0 7; 10 10 7; 0 -10 7; -10 0 7; -10 -10 7; -10 10 7; 10 -10 7;
% 0 0 8; 0 1 8; 1 0 8; 1 1 8; 0 -1 8; -1 0 8; -1 -1 8; -1 1 8; 1 -1 8;
0 0 9; 0 10 9; 10 0 9; 10 10 9; 0 -10 9; -10 0 9; -10 -10 9; -10 10 9; 10 -10 9;]'*1e-3;
voxel_targets = cart2grid(kgrid, target_pos);
%% source and sensor
source_diameter = 1e-3;
source_mag = 100;
sensor_diameter = 1e-3;
sensor_map = zeros(Nx,Ny);
for k = 1:numel(scanning_pos_x_grid)
disc_tmp = makeDisc(Nx, Ny, scanning_pos_x_grid(k), scanning_pos_y_grid(k), round(0.5*sensor_diameter/grid_interval(1)));
sensor_map = sensor_map|disc_tmp;
end
sensor_mask = zeros(Nx, Ny, Nz);
sensor_mask(:,:, Nz/2) = sensor_map;
sensor.mask = sensor_mask;
%%
medium.sound_speed = c * ones(Nx, Ny, Nz); % [m/s]
medium.sound_speed(voxel_targets==1) = 1600; % [m/s] speed of sound in hard target
medium.sound_speed(:,:,Nz/2:end) = 340; % [m/s] speed of sound in air
medium.density = 1000 * ones(Nx, Nx, Nz); % [kg/m^3] Soft tissue density
medium.density(voxel_targets==1) = 1200; % [kg/m^3] tissue density of blood
medium.density(:,:,Nz/2:end) = 1.3; % [kg/m^3] density of air
%%
% for k = 1:numel(scanning_pos_y)
for k = 1
% make source
ball_tmp = makeBall(Nx, Ny, Nz, scanning_pos_x_grid(k), scanning_pos_y_grid(k), Nz/2, round(0.5*source_diameter/grid_interval(1)));
ball_tmp(:,:,1:Nz/2-1) = 0;
source.p0 = source_mag * ball_tmp;
% run the simulation
sensorDataTmp = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
end