acousticPressure3D.m
function P = acousticPressure3D(phases,xLim,xElements,yLim,yElements,zLim,zElements)
% ACOUSTICPRESSURE3D Calculates the 3D acoustic pressure field
%  P = acousticPressure3D(phases,xLim,xElements,yLim,yElements,zLim,zElements)
%  Returns an yElements-by-xElements-by-zElements 3D acoustic pressure field
%  from 64 transducers with given phase delays. The transducers are modeled
%  as vibrating pistons.


%%% transducer constant 
P0=5; 
%%% transducer radius
r = 0.005; 
%%% emitted frequency
f = 40000; 
w = 2*pi*f; 
%%% sound velocity in the carrier fluid (air) 
c0 = 343; 
k = w/c0;

% transducers' coordinates
trans = linspace(-r*7,r*7,8);
[tx,ty] = meshgrid(trans,trans);
tx = reshape(tx,[1,1,1,64]);
ty = reshape(ty,[1,1,1,64]);

x = linspace(xLim(1),xLim(2),xElements);
y = linspace(yLim(1),yLim(2),yElements);
z = linspace(zLim(1),zLim(2),zElements);
[x,y,z] = meshgrid(x,y,z);

X = repmat(x,1,1,1,64) - repmat(tx,yElements,xElements,zElements,1);
Y = repmat(y,1,1,1,64) - repmat(ty,yElements,xElements,zElements,1);
Z = repmat(z,1,1,1,64);
phi = repmat(reshape(phases,1,1,1,64),yElements,xElements,zElements,1);

    d = sqrt(X.^2 + Y.^2 + Z.^2);
    sinTheta = sqrt(X.^2 + Y.^2)./d;
    dirFunc = 2*besselj(1, k.*r.*sinTheta)./(d.*k.*r.*sinTheta);
    dirFunc(isnan(dirFunc)) = 1;
    p = exp(1i*(d.*k+phi)).*P0.*dirFunc;
    P = abs(sum(p,4));
end