This small piece of code illustrates what I have just written (you just have to specify the correct path in the 3 first lines).

Anthony

********************************* CODE ******************************

clear all

close all

clc

infilename = 'C:\Users\Grisey\Desktop\checkSource\input.h5';

outfilename = 'C:\Users\Grisey\Desktop\checkSource\output.h5';

pathToKwave = '"C:\Users\Grisey\k-wave\k-wave-toolbox-version-1.0-cpp-windows-binaries\Windows Binaries\kspaceFirstOrder3D-OMP.exe"'; %(don't forget the ")

%% Paramètres

myPMLsize = 20;

medium.density = 1000;

medium.sound_speed = 1500;

source_mag = 1e5;

%% Grille de calcul

f = 1e5;

N = 128;

h = 1e-3;

kgrid = makeGrid(N, h, N, h, N, h);

%% Definition du transducteur

R = 10e-3;

transducerMask = zeros(N,N,N);

xgrid = h*((1:N)-N/2);

ygrid = h*((1:N)-N/2);

zgrid = h*((1:N)-N/2);

[X,Y,Z] = meshgrid(xgrid,ygrid,zgrid);

Rmap = (X.^2+Y.^2+Z.^2).^0.5;

transducerMask(Rmap<=R)=1;

%% Simulation

% create the time array

[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);

periode = 1/f;

recStartTime = numel(kgrid.t_array)-2*round(periode/dt)+1;

source.p_mask = logical(transducerMask);

source.p = source_mag*sin(2*pi*f*kgrid.t_array);

source.p_mode = 'dirichlet';

% assign the input options

input_args = {'PlotPML', true, 'PMLSize', myPMLsize,'Smooth',[false true true]};

% run the simulation

toto = ones(N, N, N);

sensor.mask = logical(toto);

clear toto;

kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:},'SaveToDisk',infilename);

status = system([pathToKwave ' -i ' infilename ' -o ' outfilename ' --p_rms --I_avg -s' num2str(recStartTime)]);

%% Lecture du résultat

p_rms = reshape(h5read(outfilename,'/p_rms'),N,N,N);

Ix_avg = reshape(h5read(outfilename,'/Ix_avg'),N,N,N);

Iy_avg = reshape(h5read(outfilename,'/Iy_avg'),N,N,N);

Iz_avg = reshape(h5read(outfilename,'/Iz_avg'),N,N,N);

H = -divergence(Ix_avg,Iy_avg,Iz_avg)/h;

% Calcul de Ir

Ianalytical = source_mag^2*R^2./(2*medium.sound_speed*medium.density.*Rmap.^2);

theta = acos(Z./Rmap);

phi = atan(Y./X)+pi*(X<0);

Ix_analytical = Ianalytical.*cos(phi).*sin(theta);

Iy_analytical = Ianalytical.*sin(phi).*sin(theta);

Iz_analytical = Ianalytical.*cos(theta);

figure(1)

subplot(1,3,1)

imagesc(squeeze(H(round(N/2),:,:)).*(1-squeeze(transducerMask(round(N/2),:,:))))

caxis([-0.1e5 2e4])

title('-div(k-wave intensity) (order 1 uncentered divergence)')

subplot(1,3,2)

Hanalytical = -divergence(Ix_analytical,Iy_analytical,Iz_analytical)/h;

imagesc(squeeze(Hanalytical(round(N/2),:,:)).*(1-squeeze(transducerMask(round(N/2),:,:))))

caxis([-0.1e5 2e4])

title('-div(analytical intensity) (order 1 uncentered divergence)')

subplot(1,3,3)

divx = (Ix_avg(3:end,2:end-1,2:end-1)-Ix_avg(1:end-2,2:end-1,2:end-1))/2/h;

divy = (Iy_avg(2:end-1,3:end,2:end-1)-Iy_avg(2:end-1,1:end-2,2:end-1))/2/h;

divz = (Iz_avg(2:end-1,2:end-1,3:end)-Iz_avg(2:end-1,2:end-1,1:end-2))/2/h;

HnumericalOrder2 = -divx-divy-divz;

imagesc(squeeze(HnumericalOrder2 (round(N/2),:,:)).*(1-squeeze(transducerMask(round(N/2),2:end-1,2:end-1))))

caxis([-0.1e5 2e4])

title('-div(k-wave intensity) (order 2 centered divergence)')