function [RP,param,SPECT,IDX] = pfield(varargin)

%PFIELD   RMS acoustic pressure field of a linear or convex array
%   RP = PFIELD(X,Y,Z,DELAYS,PARAM) returns the radiation pattern of a
%   uniform LINEAR or CONVEX array whose elements are excited at different
%   time delays (given by the vector DELAYS). The radiation pattern RP is
%   given in terms of the root-mean-square (RMS) of acoustic pressure. The
%   characteristics of the array and transmission must be given in the
%   structure PARAM (see below for details). The radiation pattern is
%   calculated at the points specified by (X,Y,Z).
%
%   RP = PFIELD(X,Z,DELAYS,PARAM) or RP = PFIELD(X,[],Z,DELAYS,PARAM)
%   disregards elevation focusing (PARAM.focus is ignored) and assumes that
%   Y=0 (2-D space). The computation is faster in 2-D.
%
%   >--- Try it: enter "pfield" in the command window for an example ---<
%
%   Units: X,Y,Z must be in m; DELAYS must be in s.
%
%   DELAYS can be a matrix. This syntax can be used to simulate MLT
%   (multi-line transmit) sequences. In this case, each ROW represents a
%   delay series. For example, to create a 4-MLT sequence with a 64-element
%   phased array, the DELAYS matrix must have 4 rows and 64 columns (size =
%   [4 64]).
%
%   PFIELD is called by SIMUS to simulate ultrasound RF radio-frequency
%   signals generated by an ultrasound uniform linear or convex array.
%
%   ---
%   NOTE #1: X-, Y-, and Z-axes
%   Conventional axes are used:
%   i)  For a LINEAR array, the X-axis is PARALLEL to the transducer and
%       points from the first (leftmost) element to the last (rightmost)
%       element (X = 0 at the CENTER of the transducer). The Z-axis is
%       PERPENDICULAR to the transducer and points downward (Z = 0 at the
%       level of the transducer, Z increases as depth increases). The
%       Y-axis is such that the coordinates are right-handed.
%   ii) For a CONVEX array, the X-axis is parallel to the chord and Z = 0
%       at the level of the chord.
%   ---
%   NOTE #2: Simplified method: Directivity
%   By default, the calculation is made faster by assuming that the
%   directivity of the elements is dependent only on the central frequency.
%   This simplification very little affects the pressure field in most
%   situations (except near the array). To turn off this option, use
%   OPTIONS.FullFrequencyDirectivity = true.
%   (see ADVANCED OPTIONS below).
%   ---
%
%   PARAM is a structure that contains the following fields:
%   -------------------------------------------------------
%       *** TRANSDUCER PROPERTIES ***
%   1)  PARAM.fc: central frequency (in Hz, REQUIRED)
%   2)  PARAM.pitch: pitch of the array (in m, REQUIRED)
%   3)  PARAM.width: element width (in m, REQUIRED)
%        or PARAM.kerf: kerf width (in m, REQUIRED)
%        note: width = pitch-kerf
%   4)  PARAM.focus: elevation focus (in m, ignored if Y is not given)
%            The default is Inf (no elevation focusing)
%   5)  PARAM.height: element height (in m, ignored if Y is not given)
%            The default is Inf (no elevation focusing)
%   6)  PARAM.radius: radius of curvature of the probe (in m)
%            The default is Inf (rectilinear array)
%   7)  PARAM.bandwidth: pulse-echo 6dB fractional bandwidth (in %)
%            The default is 75%.
%   8)  PARAM.baffle: property of the baffle:
%            'soft' (default), 'rigid', or a scalar > 0.
%            See "Note on BAFFLE properties" below for details
%
%       *** MEDIUM PARAMETERS ***
%   9)  PARAM.c: longitudinal velocity (in m/s, default = 1540 m/s)
%   10) PARAM.attenuation: attenuation coefficient (dB/cm/MHz, default: 0)
%            Notes: A linear frequency-dependence is assumed.
%                   A typical value for soft tissues is ~0.5 dB/cm/MHz.
%
%       *** TRANSMIT PARAMETERS ***
%   11) PARAM.TXapodization: transmit apodization (default: no apodization)
%   12) PARAM.TXnow: number of wavelengths of the TX pulse (default: 1)
%   13) PARAM.TXfreqsweep: frequency sweep for a linear chirp (default: [])
%                          To be used to simulate a linear TX chirp.
%
%   Other syntaxes:
%   --------------
%   i}  [RP,PARAM] = PFIELD(...) also returns the complete list of
%       parameters including the default values.
%   ii} [...] = PFIELD without any input argument provides an interactive
%       example designed to produce a focused ultrasound beam using a 2.7
%       MHz phased-array transducer (without elevation focusing).
%
%
%   Note on CHIRP signals:
%   ---------------------
%   Linear chirps are characterized by PARAM.TXnow, PARAM.fc and
%   PARAM.TXfreqsweep. The transmitted pulse The transmitted pulse has a
%   duration of approximately T (= PARAM.TXnow/PARAM.fc), with the
%   amplitude and phase defined over the time interval -T/2 to +T/2. The
%   total frequency sweep is DeltaF (= PARAM.TXfreqsweep): the frequencies
%   changes linearly from (PARAM.fc - DeltaF/2) to (PARAM.fc + DeltaF/2)
%   in the defined time interval.
%
%
%   Note on BAFFLE property:
%   -----------------------
%   In PFIELD, it is assumed by default that the array elements are
%   embedded in an infinite SOFT baffle. To modify the property of the
%   baffle, modify the field PARAM.baffle:
%       1) 'rigid'
%       2) 'soft' (this is the default)
%       3) a nonnegative scalar Alpha,
%          with Alpha = (medium impedance)/(baffle impedance)
%          Note: Alpha = 0 => 'rigid'; Alpha >> 1 => 'soft'
%
%   The baffle property affects the obliquity factor included in the
%   directivity of the elements. This obliquity factor is not 1 if the
%   baffle is not rigid. A general case (see case #3 below) can be chosen
%   by specifying an impedance ratio. For details, refer to the
%   corresponding papers.
%   1) For a rigid baffle => obliquity factor = 1.
%   2) For a soft baffle => obliquity factor = cos(Theta).
%      Selfridge et al. Appl Phys Lett 37(1), 35-36 (1980)
%      "A theory for the radiation pattern of a narrow-strip acoustic
%      transducer." <a
%      href="matlab:web('http://scitation.aip.org/content/aip/journal/apl/37/1/10.1063/1.91692')">Paper here</a>
%   3) General baffle => obliquity factor = cos(Theta)/(cos(Theta)+Alpha)
%      with Alpha = (medium impedance)/(baffle impedance).
%      Pesqué et al. IEEE Ultrasonics Symposium, (1984)
%      "Effect of the planar baffle impedance in acoustic radiation of a
%      phased array element theory and experimentation." <a
%      href="matlab:web('http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=1535402')">Paper here</a>
%      Example: For a baffle of impedance 2.8 MRayl (epoxy) adjacent to
%               soft tissues of impedance 1.6 MRayls, alpha = 1.75.
%
%
%   ADVANCED OPTIONS:
%   ----------------
%       %-- FREQUENCY SAMPLES --%
%   1)  Only frequency components of the transmitted signal in the range
%       [0,2fc] with significant amplitude are considered. The default
%       relative amplitude is -60 dB in PFIELD. You can change this value
%       by using the following:
%           [...] = PFIELD(...,OPTIONS),
%       where OPTIONS.dBThresh is the threshold in dB (default = -60).
%   ---
%       %-- FULL-FREQUENCY DIRECTIVITY --%   
%   2)  By default, the directivity of the elements depends on the center
%       frequency only. This makes the algorithm faster. To make the
%       directivities fully frequency-dependent, use: 
%           [...] = PFIELD(...,OPTIONS),
%       with OPTIONS.FullFrequencyDirectivity = true (default = false).
%   ---
%       %-- ELEMENT SPLITTING --%   
%   3)  Each transducer element of the array is split into small segments.
%       The length of these small segments must be small enough to ensure
%       that the far-field model is accurate. By default, the elements are
%       split into M segments, with M being defined by:
%           M = ceil(element_width/smallest_wavelength);
%       To modify the number M of subelements by splitting, you may adjust
%       OPTIONS.ElementSplitting. For example, OPTIONS.ElementSplitting = 1
%   ---
%       %-- WAIT BAR --%   
%   4)  If OPTIONS.WaitBar is true, a wait bar appears (only if the number
%       of frequency samples >10). Default is true.
%   ---
%
%
%   Notes regarding the model & REFERENCES:
%   --------------------------------------
%   1) PFIELD works for uniform linear or convex arrays. A uniform array
%      has identical elements along a rectilinear or curved line in space
%      with uniform spacing. Each element is split into small segments (if
%      required). The radiation patterns in the x-z plane are derived by
%      using Fraunhofer (far-field) equations. Those in the x-y elevational
%      plane are derived by using a Fresnel (paraxial) approximation.
%   2) The paper that describes the first 2-D version of PFIELD is:
%      SHAHRIARI S, GARCIA D. Meshfree simulations of ultrasound vector
%      flow imaging using smoothed particle hydrodynamics. Phys Med Biol,
%      2018;63:205011. <a
%      href="matlab:web('https://www.biomecardio.com/publis/physmedbio18.pdf')">PDF here</a>
%   3) The paper that describes the theory of the full (2-D + 3-D) version
%      of PFIELD is under review.
%
%
%   EXAMPLES:
%   --------
%
%   %-- Example #1: 2-D focused pressure field with a cardiac phased-array
%   % Phased-array @ 2.7 MHz:
%   param = getparam('P4-2v');
%   % Focus position:
%   x0 = 2e-2; z0 = 5e-2;
%   % TX time delays:
%   dels = txdelay(x0,z0,param);
%   % Grid:
%   x = linspace(-4e-2,4e-2,200);
%   z = linspace(0,10e-2,200);
%   [x,z] = meshgrid(x,z);
%   % RMS pressure field:
%   P = pfield(x,z,dels,param);
%   imagesc(x(1,:),z(:,1),20*log10(P/max(P(:))))
%   colormap hot, axis equal tight
%   hold on, plot(x0,z0,'k*'), hold off
%   title('2-D focused field ([-20,0] dB)')
%   caxis([-20 0])
%
%   %-- Example #2: 3-D focused pressure field with a cardiac phased-array
%   % Phased-array @ 2.7 MHz:
%   param = getparam('P4-2v');
%   % Focus position:
%   x0 = 0; z0 = 5e-2;
%   % TX time delays:
%   dels = txdelay(x0,z0,param);
%   % x-, y-, z-coordinates:
%   x = linspace(-4e-2,4e-2,200);
%   y = linspace(-0.75,.75,50)*param.height;
%   z = linspace(0,10e-2,200);
%   % RMS pressure field, azimuthal plane:
%   [xaz,zaz] = meshgrid(x,z);
%   yaz = zeros(size(xaz));
%   Paz = pfield(xaz,yaz,zaz,dels,param);
%   % RMS pressure field, elevation plane:
%   [yel,zel] = meshgrid(y,z);
%   xel = zeros(size(yel));
%   Pel = pfield(xel,yel,zel,dels,param);
%   % RMS pressure field, focal plane:
%   [xfo,yfo] = meshgrid(x,y);
%   zfo = ones(size(xfo))*param.focus;
%   Pfo = pfield(xfo,yfo,zfo,dels,param);
%   % Figure:
%   Pmax = max(Paz(:));
%   surf(xaz,yaz,zaz,20*log10(Paz/Pmax)), shading flat
%   hold on
%   surf(xel,yel,zel,20*log10(Pel/Pmax)), shading flat
%   surf(xfo,yfo,zfo,20*log10(Pfo/Pmax)), shading flat
%   hold off
%   % Fine-tune the figure:
%   axis equal
%   set(gca,'ZDir','reverse')
%   title('3-D focused field ([-20,0] dB)')
%   caxis([-20 0])
%   alpha color
%
%   %-- Example #3: 3-D MLT multi-line transmission
%   % Phased-array @ 2.7 MHz:
%   param = getparam('P4-2v');
%   % Delays for a MLT multi-line transmission:
%   x0 = 2e-2; z0 = 5e-2;
%   dels = txdelay([-x0 0 x0],[z0 sqrt(x0^2+z0^2) z0],param);
%   % Grid:
%   x = linspace(-4e-2,4e-2,200);
%   z = linspace(0,10e-2,200);
%   [x,z] = meshgrid(x,z);
%   % RMS pressure field:
%   P = pfield(x,zeros(size(x)),z,dels,param);
%   imagesc(x(1,:),z(:,1),20*log10(P/max(P(:))))
%   colormap hot, axis equal tight
%   title('3-D MLT field ([-20,0] dB)')
%   caxis([-20 0])
%
%
%   This function is part of <a
%   href="matlab:web('https://www.biomecardio.com/MUST')">MUST</a> (Matlab UltraSound Toolbox).
%   MUST (c) 2020 Damien Garcia, LGPL-3.0-or-later
%
%   See also TXDELAY, SIMUS, GETPARAM, GETPULSE.
%
%   -- Damien Garcia -- 2015/02, last update 2022/04/07
%   website: <a
%   href="matlab:web('https://www.biomecardio.com')">www.BiomeCardio.com</a>


MyFolder = 'C:\Users\garcia\Documents\-- Sync Home'; % PLZ do not modify,
% unless you are getting tired of the advertising message...
% In this case, use any folder that you own.
% 
if ~exist(MyFolder,'dir')
    mlock % NOTE: use MUNLOCK if you need to clear PFIELD from memory
    persistent pfieldLastUse %#ok
end


if nargin==0
    if nargout>0
        [RP,param] = RunTheExample;
    else
        RunTheExample;
    end
    return
end
narginchk(4,7)

%-- Check if we have the 'quick' option
% masked option: PFIELD(...,'quick') for quick tests
if strcmpi(varargin{end},'quick')
    NArg = nargin-1;
    isQUICK = true;
else
    NArg = nargin;
    isQUICK = false;
end

%-- Input variables: X,Y,Z,DELAYS,PARAM,OPTIONS
x = varargin{1};
switch NArg
    case 4 % pfield(X,Z,DELAYS,PARAM)
        y = [];
        z = varargin{2};
        delaysTX = varargin{3};
        param = varargin{4};
        assert(isstruct(param),'PARAM must be a structure.')
        options = [];
    case 5
        if isstruct(varargin{4}) % pfield(X,Z,DELAYS,PARAM,OPTIONS)
            y = [];
            z = varargin{2};
            delaysTX = varargin{3};
            param = varargin{4};
            options = varargin{5};
        else % pfield(X,Y,Z,DELAYS,PARAM)
            y = varargin{2};
            z = varargin{3};
            delaysTX = varargin{4};
            param = varargin{5};
            options = [];
        end
    otherwise % pfield(X,Y,Z,DELAYS,PARAM,OPTIONS)
            y = varargin{2};
            z = varargin{3};
            delaysTX = varargin{4};
            param = varargin{5};
            options = varargin{6};
end

%-- 'quick' options
if isQUICK
    options.dBThresh = -20;
    options.ElementSplitting = 1;
    options.FullFrequencyDirectivity = false;
    options.FrequencyStep = 1.5;
end

%-- Elevation focusing and X,Y,Z size
if isempty(y)
    ElevationFocusing = false;
    assert(isequal(size(x),size(z)),...
        'X and Z must be of same size.')
    y = 0;
else
    ElevationFocusing = true;
    assert(isequal(size(x),size(y),size(z)),...
        'X, Y, and Z must be of same size.')
end

%-- Check the transmit delays
assert(isnumeric(delaysTX) && all(delaysTX(~isnan(delaysTX))>=0),...
    'DELAYS must be a nonnegative array.')
if isvector(delaysTX) % we need a row vector
    delaysTX = delaysTX(:).'; 
end
NumberOfElements = size(delaysTX,2);
% Note: param.Nelements can be required in other functions of the
%       Matlab Ultrasound Toolbox
if isfield(param,'Nelements')
    assert(param.Nelements==NumberOfElements,...
        ['DELAYS must be of length PARAM.Nelements.',13,...
        'If DELAYS is a matrix, its number of rows must be PARAM.Nelements.'])
elseif nargout>1
    param.Nelements = NumberOfElements;
end
%--
% Note: delaysTX can be a matrix. This option can be used for MLT
% (multi-line transmit) for example. In this case, each row represents a
% delay series. For example, for a 4-MLT sequence with a 64-element phased
% array, delaysTX has 4 rows and 64 columns, i.e. size(delaysTX) = [4 64].
%--
delaysTX = delaysTX.';

%-- Check if PFIELD is called by SIMUS or MKMOVIE
isSIMUS = false;
isMKMOVIE = false;
if isfield(options,'CallFun')
    isSIMUS = strcmpi(options.CallFun,'simus');
    isMKMOVIE = strcmpi(options.CallFun,'mkmovie');
end

%-- Check if SIMUS is running on a parallel pool
try
    onppool = isSIMUS && options.ParPool;
catch
    onppool = false;
end

%-- Advertising message & Statistics
if ~onppool && ~exist(MyFolder,'dir')
    if isempty(pfieldLastUse) || (now-pfieldLastUse)>1
        pfieldLastUse = now;
        
        % An ad appears once a day (once PFIELD is used) in a small message
        % box. Please do not remove it!
        AdMessage
        
        % -- Statistics on daily use -- The number of uses is incremented
        % daily, and the information is sent in an FTP file for statistical
        % purpose. Don't worry, there is no spam! The MUSTstat file is
        % encrypted (p code) because it contains a password. Please do not
        % remove it!
        try MUSTstat, catch, end
    end
end


%---------------------------%
% Check the PARAM structure %
%---------------------------%

param = IgnoreCaseInFieldNames(param);

%-- 1) Center frequency (in Hz)
assert(isfield(param,'fc'),...
    'A center frequency value (PARAM.fc) is required.')
fc = param.fc; % central frequency (Hz)

%-- 2) Pitch (in m)
assert(isfield(param,'pitch'),'A pitch value (PARAM.pitch) is required.')
pitch = param.pitch;

%-- 3) Element width and/or Kerf width (in m)
if isfield(param,'width') && isfield(param,'kerf')
    assert(abs(pitch-param.width-param.kerf)<eps,...
        'The pitch must be equal to (kerf width + element width).')
elseif isfield(param,'kerf')
    param.width = pitch-param.kerf;
elseif isfield(param,'width')
    param.kerf = pitch-param.width;
else
    error(['An element width (PARAM.width) ',...
        'or kerf width (PARAM.kerf) is required.'])
end
ElementWidth = param.width;

%-- 4) Elevation focus (in m)
if ~isfield(param,'focus')
    param.focus = Inf; % default = no elevation focusing
end
Rf = param.focus;
assert(isnumeric(Rf) && isscalar(Rf) && Rf>0,...
    'The element focus must be positive.') 

%-- 5) Element height (in m)
if ~isfield(param,'height')
    param.height = Inf; % default = line array
end
ElementHeight = param.height;
assert(isnumeric(ElementHeight) && isscalar(ElementHeight) &&...
    ElementHeight>0,'The element height must be positive.')

%-- 6) Radius of curvature (in m) - convex array
if ~isfield(param,'radius')
    param.radius = Inf; % default = linear array
end
RadiusOfCurvature = param.radius;
assert(isnumeric(RadiusOfCurvature) && isscalar(RadiusOfCurvature) &&...
    RadiusOfCurvature>0,'The radius of curvature must be positive.')

%-- 7) Fractional bandwidth at -6dB (in %)
if ~isfield(param,'bandwidth')
    param.bandwidth = 75;
end
assert(param.bandwidth>0 && param.bandwidth<200,...
    'The fractional bandwidth at -6 dB (PARAM.bandwidth, in %) must be in ]0,200[')

%-- 8) Baffle
%   An obliquity factor will be used if the baffle is not rigid
%   (default = SOFT baffle)
if ~isfield(param,'baffle')
    param.baffle = 'soft'; % default
end
if strcmpi(param.baffle,'rigid')
    NonRigidBaffle = false;
elseif strcmpi(param.baffle,'soft')
    NonRigidBaffle = true;
elseif isscalar(param.baffle)
    assert(param.baffle>0,...
        'The ''baffle'' field scalar must be positive')
    NonRigidBaffle = true;
else
    error('The ''baffle'' field must be ''rigid'',''soft'' or a positive scalar')
end

%-- 9) Longitudinal velocity (in m/s)
if ~isfield(param,'c')
    param.c = 1540; % default value
end
c = param.c; % speed of sound (m/s)

%-- 10) Attenuation coefficient (in dB/cm/MHz)
if ~isfield(param,'attenuation') % no attenuation, alpha_dB = 0
    param.attenuation = 0;
    alpha_dB = 0;
else
    alpha_dB = param.attenuation;
    assert(isscalar(alpha_dB) && isnumeric(alpha_dB) && alpha_dB>=0,...
        'PARAM.attenuation must be a nonnegative scalar')
end

%-- 11) Transmit apodization (no unit)
if ~isfield(param,'TXapodization')
    param.TXapodization = ones(1,NumberOfElements);
else
    assert(isvector(param.TXapodization) && isnumeric(param.TXapodization),...
        'PARAM.TXapodization must be a vector')
    assert(numel(param.TXapodization)==NumberOfElements,...
        'PARAM.TXapodization must be of length = (number of elements)')
end
% apodization is 0 where TX delays are NaN:
idx = isnan(delaysTX);
param.TXapodization(any(idx,2)) = 0;
delaysTX(idx) = 0;

%-- 12) TX pulse: Number of wavelengths
if ~isfield(param,'TXnow')
    param.TXnow = 1;
end
NoW = param.TXnow;
assert(isscalar(NoW) && isnumeric(NoW) && NoW>0,...
    'PARAM.TXnow must be a positive scalar.')

%-- 13) TX pulse: Frequency sweep for a linear chirp
if ~isfield(param,'TXfreqsweep') || isinf(NoW)
    param.TXfreqsweep = [];
end
FreqSweep = param.TXfreqsweep;
assert(isempty(FreqSweep) ||...
    (isscalar(FreqSweep) && isnumeric(FreqSweep) && FreqSweep>0),...
    'PARAM.TXfreqsweep must be empty (windowed sine) or a positive scalar (linear chirp).')

%----------------------------------%
% END of Check the PARAM structure %
%----------------------------------%



%-----------------------------%
% Check the OPTIONS structure %
%-----------------------------%

options = IgnoreCaseInFieldNames(options);

%-- 1) dB threshold
%     (in dB: faster computation if lower value, but less accurate)
if ~isfield(options,'dBThresh')
    options.dBThresh = -60; % default is -60dB in PFIELD
end
assert(isscalar(options.dBThresh) && isnumeric(options.dBThresh) &&...
    options.dBThresh<=0,'OPTIONS.dBThresh must be a nonpositive scalar.')

%-- 2) Frequency-dependent directivity?
if isfield(options,'FullFrequencyDirectivity')
    isFFD = options.FullFrequencyDirectivity;
else
    isFFD = false; % default
    % By default, the directivity of the elements depends on the center
    % frequency only. This makes the algorithm faster. 
end
assert(isscalar(isFFD) && islogical(isFFD),...
    'OPTIONS.FullFrequencyDirectivity must be a logical scalar (true or false).')

%-- 3) Element splitting
%
% --- A short note about the algorithm:
% A far-field model is used in PFIELD. Each transducer element of the array
% is split into M small segments, so that these M small segments are
% smaller than one wavelength. The far-field condition is acceptable for
% these small segments.
%---
if isfield(options,'ElementSplitting') && ~isempty(options.ElementSplitting)
    M = options.ElementSplitting;
    assert(isscalar(M) & M==round(M) & M>0,...
        'OPTIONS.ElementSplitting must be a positive integer.')
else
    LambdaMin = c/(fc*(1+param.bandwidth/200));
    M = ceil(ElementWidth/LambdaMin);
end

%-- 4) Wait bar
if ~isfield(options,'WaitBar')
    options.WaitBar = true;
end
assert(isscalar(options.WaitBar) && islogical(options.WaitBar),...
    'OPTIONS.WaitBar must be a logical scalar (true or false).')

%-- Advanced (masked) options: Frequency step (scaling factor)
% The frequency step is determined automatically. It is tuned to avoid
% significant interferences due to unadapted discretization. The frequency
% step can also be adjusted by using a scaling factor. For a fast check,
% you may use a scaling factor>1. For a smoother result, you may use a
% scaling factor<1.
if ~isfield(options,'FrequencyStep')
    options.FrequencyStep = 1;
end
assert(isscalar(options.FrequencyStep) &&...
    isnumeric(options.FrequencyStep) && options.FrequencyStep>0,...
    'OPTIONS.FrequencyStep must be a positive scalar.')

%------------------------------------%
% END of Check the OPTIONS structure %
%------------------------------------%


%-----
% SIMUS and MKMOVIE first run a PFIELD with empty X,Y,Z to detect possible
% syntax errors.
if isSIMUS && isempty(x), RP = []; return, end
if isMKMOVIE && isempty(x) && nargout==2, RP = []; SPECT = []; return, end
%-----


%------------------------------------%
% POINT LOCATIONS, DISTANCES & GRIDS %
%------------------------------------%

siz0 = size(x);
nx = numel(x);

%-- Coordinates of the points where pressure is needed
x = x(:); y = y(:); z = z(:);
if isMKMOVIE
    x = [x; options.x(:)];
    z = [z; options.z(:)];
    % Note with MKMOVIE:
    % We must consider the points of the image grid + the points of the
    % scatterers (if any). The scatterer coordinates are in options.x and
    % options.z.
    % Note: there is no elevation focusing with MKMOVIE (2-D only).
end
% cast x, y, and z to single class
x = cast(x,'single');
y = cast(y,'single');
z = cast(z,'single');

%-- Centers of the tranducer elements (x- and z-coordinates)
if isinf(RadiusOfCurvature)
    % we have a LINEAR ARRAY
    xe = ((0:NumberOfElements-1)-(NumberOfElements-1)/2)*param.pitch;
    ze = zeros(1,NumberOfElements);
    THe = zeros(1,NumberOfElements);
else
    % we have a CONVEX ARRAY
    chord = 2*RadiusOfCurvature*...
        sin(asin(param.pitch/2/RadiusOfCurvature)*(NumberOfElements-1));
    h = sqrt(RadiusOfCurvature^2-chord^2/4); % apothem
    % https://en.wikipedia.org/wiki/Circular_segment
    % THe = angle of the normal to element #e with respect to the z-axis
    THe = linspace(atan2(-chord/2,h),atan2(chord/2,h),NumberOfElements);
    ze = RadiusOfCurvature*cos(THe);
    xe = RadiusOfCurvature*sin(THe);
    ze = ze-h;
end


%-- Centroids of the sub-elements
%-- note: The elements have been split into M sub-elements.
% X-position (xi) and Z-position (zi) of the centroids of the sub-elements
% (relative to the centers of the transducer elements).
% The values in xi,zi are in the range ]-ElementWidth/2 ElementWidth/2[
% (if M=1, then xi = zi = 0 for a rectilinear array).
SegLength = ElementWidth/M;
tmp = -ElementWidth/2 + SegLength/2 + (0:M-1)*SegLength;
xi = reshape(tmp,[1 1 M]).*cos(THe);
zi = reshape(tmp,[1 1 M]).*sin(-THe);

%-- Out-of-field points
% Null pressure will be assigned to out-of-field points.
isOUT = z<0;
if isfinite(RadiusOfCurvature)
    isOUT = isOUT |...
        (x.^2+(z+h).^2)<=RadiusOfCurvature^2;
end

%-- Variables that we need:
%
% Note: r = distance between the segment centroid and the point of interest
%       Th = angle relative to the element normal axis.
%       sinT = sine of the angle relative to the element normal axis.
%       r, Th, and sinT are of size [numel(x) NumberOfElements M].
%
dxi = x-xi-xe;
d2 = dxi.^2+(z-zi-ze).^2;
%---
% we'll have 1/sqrt(r) or 1/r:
% small d2 values are replaced by (lambda/4)^2
smallD2 = (c/fc/4)^2;
d2(d2<smallD2) = smallD2;
%---
r = sqrt(d2+y.^2);
Th = asin(dxi./sqrt(d2))-THe;
sinT = sin(Th);
clear dxi d2

%-------------------------------------------%
% end of POINT LOCATIONS, DISTANCES & GRIDS %
%-------------------------------------------%


mysinc = @(x) sin(abs(x)+eps)./(abs(x)+eps); % cardinal sine
% [note: In MATLAB, sinc is sin(pi*x)/(pi*x)]

%-------------------%
% FREQUENCY SPECTRA %
%-------------------%

%-- FREQUENCY SPECTRUM of the transmitted pulse
if isempty(FreqSweep)
    % We want a windowed sine of width PARAM.TXnow
    T = NoW/fc; % temporal pulse width
    if isinf(T); T = 1e6; end
    wc = 2*pi*fc;
    pulseSpectrum = @(w) 1i*(mysinc(T*(w-wc)/2)-mysinc(T*(w+wc)/2));
else
    % We want a linear chirp of width PARAM.TXnow
    % (https://en.wikipedia.org/wiki/Chirp_spectrum#Linear_chirp)
    T = NoW/fc; % temporal pulse width
    if isinf(T); T = 1e6; end
    wc = 2*pi*fc;
    dw = 2*pi*FreqSweep;
    s2 = @(w) sqrt(pi*T/dw)*exp(-1i*(w-wc).^2*T/2/dw).*...
        (fresnelint((dw/2+w-wc)/sqrt(pi*dw/T)) + fresnelint((dw/2-w+wc)/sqrt(pi*dw/T)));
    pulseSpectrum = @(w) (1i*s2(w)-1i*s2(-w))/T;
end

%-- FREQUENCY RESPONSE of the ensemble PZT + probe
% We want a generalized normal window (6dB-bandwidth = PARAM.bandwidth)
% (https://en.wikipedia.org/wiki/Window_function#Generalized_normal_window)
wB = param.bandwidth*wc/100; % angular frequency bandwidth
p = log(126)/log(2*wc/wB); % p adjusts the shape
probeSpectrum = @(w) exp(-(abs(w-wc)/(wB/2/log(2)^(1/p))).^p);
% The frequency response is a pulse-echo (transmit + receive) response. A
% square root is thus required when calculating the pressure field:
probeSpectrum = @(w) sqrt(probeSpectrum(w));
% Note: The spectrum of the pulse (pulseSpectrum) will be then multiplied
% by the frequency-domain tapering window of the transducer (probeSpectrum)

%-- FREQUENCY STEP
if isSIMUS || isMKMOVIE % PFIELD has been called by SIMUS or MKMOVIE
    df = options.FrequencyStep;
else % We are in PFIELD only (i.e. not called by SIMUS or MKMOVIE)
    % The frequency step df is chosen to avoid interferences due to
    % inadequate discretization.
    % -- df = frequency step (must be sufficiently small):
    % One has exp[-i(k r + w delay)] = exp[-2i pi(f r/c + f delay)] in the Eq.
    % One wants: the phase increment 2pi(df r/c + df delay) be < 2pi.
    % Therefore: df < 1/(r/c + delay).
    df = 1/(max(r(:)/c) + max(delaysTX(:)));
    df = options.FrequencyStep*df;
    % note: df is here an upper bound; it will be recalculated below
end

%-- FREQUENCY SAMPLES
Nf = 2*ceil(param.fc/df)+1; % number of frequency samples
f = linspace(0,2*param.fc,Nf); % frequency samples
df = f(2); % update the frequency step
%- we keep the significant components only by using options.dBThresh
S = abs(pulseSpectrum(2*pi*f).*probeSpectrum(2*pi*f));
GdB = @(f) 20*log10(S/max(S)); % gain in dB
IDX = GdB(f)>options.dBThresh;
IDX = IDX(find(IDX,1):find(IDX,1,'last'));
f = f(IDX);
nSampling = length(f);

%-- For MKMOVIE only: IDX is required in MKMOVIE
if isMKMOVIE && isempty(x), RP = []; SPECT = []; return, end

%-- we need VECTORS
pulseSPECT = pulseSpectrum(2*pi*f); % pulse spectrum
probeSPECT = probeSpectrum(2*pi*f); % probe response

%--------------------------%
% end of FREQUENCY SPECTRA %
%--------------------------%



%-- Wait bar
options.WaitBar = options.WaitBar & (nSampling>10);
if options.WaitBar
    if isSIMUS
        wbtitle = 'Let SIMUS do the work for you...';
        wbname = 'SIMUS / www.biomecardio.com';
    elseif isMKMOVIE
        wbtitle = 'Let MKMOVIE do the work for you...';
        wbname = 'MKMOVIE / www.biomecardio.com';
    else
        wbtitle = 'Let PFIELD do the work for you...';
        wbname = 'PFIELD / www.biomecardio.com';
    end
    tstart = tic;
    hwb = waitbar(0,wbtitle,'Name',wbname);
end

%-- Initialization
RP = 0; % RP = Radiation Pattern
if isSIMUS
    %- For SIMUS only (we need the full spectrum of RX signals):
    SPECT = zeros([nSampling NumberOfElements siz0(2)],'like',single(1i));
elseif isMKMOVIE
    %- For MKMOVIE only (we need the full spectrum of the pressure field):
    SPECT = zeros([nSampling nx],'like',single(1i));
elseif nargout==3
    SPECT = zeros([nSampling nx],'like',single(1i));
end

%-- Obliquity factor (baffle property)
%   An obliquity factor is required if the baffle is not rigid.
%   [Th = angle relative to the element normal axis]
if NonRigidBaffle
    if strcmpi(param.baffle,'soft')
        ObliFac = cos(Th);
    else % param.baffle is a scalar
        ObliFac = cos(Th)./(cos(Th)+param.baffle);
    end
else % 1 if rigid baffle
    ObliFac = ones(size(Th));
end
ObliFac(abs(Th)>=pi/2) = eps('single');
clear Th

%-- Note on Attenuation
% Reference: Diagnostic ultrasound imaging - inside out (T.L. Szabo)
%            Chapter 4: Attenuation
% Key reference: Acoustics for ultrasound imaging (Ben Cox, 2013)
%                Chapter 5: Acoustic attenuation and absorption
% We will use this attenuation-based wavenumber:
%   kwa = alpha_dB/8.69*f(k)/1e6*1e2; % P(z,f) = P0 exp(-alpha*f*z/8.69)
%   note: 20/log(10) ~ 8.69

%-- EXPONENTIAL arrays of size [numel(x) NumberOfElements M]
kw = 2*pi*f(1)/c; % wavenumber
kwa = alpha_dB/8.69*f(1)/1e6*1e2; % attenuation-based wavenumber
EXP = exp(-kwa*r + 1i*mod(kw*r,2*pi)); % faster than exp(-kwa*r+1i*kw*r)
%-- Exponential array for the increment wavenumber dk
dkw = 2*pi*df/c;
dkwa = alpha_dB/8.69*df/1e6*1e2;
EXPdf = exp((-dkwa + 1i*dkw)*r);

%-- We replace EXP by EXP.*ObliFac./r or EXP.*ObliFac./sqrt(r)
if ElevationFocusing
    EXP = EXP.*ObliFac./r;
    rm = mean(r,3); 
else
    EXP = EXP.*ObliFac./sqrt(r);
end
clear ObliFac r

%-- TX apodization
APOD = param.TXapodization(:);

%-- If PFIELD is called by MKMOVIE
%   Other matrices are needed if scatterers are present. If scatterers are
%   added in the movie, it is needed to simulate the propagation of the
%   backscattered echoes. We thus need the distances between the scatterers
%   and grid-points, and the corresponding EXP_RC matrix.
if isMKMOVIE && ~isempty(options.RC)
    dx = x(1:nx).'-options.x(:);
    dz = z(1:nx).'-options.z(:);
    r_RC = sqrt(dx.^2+dz.^2);
    clear dx dz
    
    % EXP_RC = exp((-kwa+1i*kw)*r_RC);
    EXP_RC = exp(-kwa*r_RC + 1i*mod(kw*r_RC,2*pi));
    EXPdf_RC = exp((-dkwa + 1i*dkw)*r_RC);
end

%-- Simplified directivity (if not dependent on frequency)
% In the "simplified directivity" version, the directivity of the elements
% depends on the center frequency ONLY. It is thus not needed to calculate
% the directivity array (DIR) in the following for-loop. This directivity
% DIR is included in the variable EXP to reduce storage.
if ~isFFD
    kc = 2*pi*fc/c; % center wavenumber
    DIR = mysinc(kc*SegLength/2*sinT); % directivity of each segment
    EXP = EXP.*DIR;
    clear DIR
end

%-- Multi-Gaussian beam model (MGBM)
% A multi-Gaussian beam model is used when integrating the elevation
% focusing.
if ElevationFocusing
    [A,B] = MGBMcoeff(4);
    % Note: The rectangular function is decomposed into a number of
    %       Gaussian functions, defined by A(n)*exp(-B(n)*y^2), with A(n)
    %       and B(n) being complex coefficients. Four coefficients offers
    %       a good compromise.
    alpha = 1i/2*(1/Rf-1./rm);
    gamma = 1i*y.^2/2./rm;
    beta2 = complex(-(y./rm).^2,0);
    clear rm
    Nmgbm = max(3,round(nSampling/20));
    k4mgbm = round(linspace(1,nSampling,Nmgbm));
    % (the MGBM field won't be calculated at each step)
end


%-----------------------------%
% SUMMATION OVER THE SPECTRUM %
%-----------------------------%

for k = 1:nSampling

    kw = 2*pi*f(k)/c; % wavenumber
    
    %--- MGBM = multi-Gaussian beam model
    % The MGBM is used to simulate the focused beams in the elevation
    % direction. These lines are run if ElevationFocusing is true.
    if ElevationFocusing && ismember(k,k4mgbm)
        % (the MGBM field is not recalculated at each step)
        MGBM = 0;
        if all(abs(y)<eps('single'))
            for n = 1:numel(A)
                MGBM = MGBM + A(n)*...
                    sqrt(pi./(kw*alpha + B(n)/ElementHeight^2));
            end
        else
            for n = 1:numel(A)
                tmp = 1./(kw*alpha + B(n)/ElementHeight^2);
                MGBM = MGBM + A(n)*sqrt(pi*tmp).*...
                    exp(kw^2*beta2/4.*tmp + kw*gamma);
            end
        end
    end
    %{
    if ElevationFocusing
        P = kw/2*(rm-Rf)./(rm*Rf+eps); % instead of: P = kw/2*(1/Rf-1./rm);
        Q = kw*y./rm;
        R = kw*y.^2/2./rm;
        MGBM = (-1)^.25*sqrt(pi)*exp(1i*(4*P.*R+Q.^2)/4./P)/2./sqrt(P).*...
            (myerfi((-1)^(3/4)*(Q-ElementHeight*P)/2./sqrt(P),64) -...
            myerfi((-1)^(3/4)*(Q+ElementHeight*P)/2./sqrt(P),64));
        whos MGBM
    end
    %}
    
    %-- Exponential array of size [numel(x) NumberOfElements M]
    % For all k, we need: EXP = exp((-kwa+1i*kw)*r)
    %                         with kw = 2*pi*f(k)/c;
    %                     and with kwa = alpha_dB/8.7*f(k)/1e6*1e2;
    % Since f(k) = f(1)+(k-1)df, we use the following recursive product:
    if k>1
        EXP = EXP.*EXPdf;
        % If PFIELD is called by MKMOVIE and if scatterers are present:
        if isMKMOVIE&&~isempty(options.RC), EXP_RC = EXP_RC.*EXPdf_RC; end
    end
        
    %-- Directivity (if frequency-dependent)
    if isFFD % isFFD = true -> frequency-dependent directivity
        DIR = kw*SegLength/2*sinT;
        DIR = mysinc(DIR);
    end
        
    %-- Radiation patterns of the single elements
    % They are the combination of the far-field patterns of the M small
    % segments that make up the single elements
    %--
    if isFFD % isFFD = true -> frequency-dependent directivity
        RPmono = mean(DIR.*EXP,3); % summation over the M small segments
    else % isFFD = false: the directivity depends on center frequency only
         % note: the directivity (DIR) has already been included in EXP
        if M>1            
            RPmono = mean(EXP,3); % summation over the M small segments
        else
            RPmono = EXP;
        end
    end
    
    %-- Complete the radiation patterns of the single elements by including
    %   the elevation focusing, if required
    if ElevationFocusing
        RPmono = RPmono.*MGBM;
    end
        
    %-- Transmit delays + Transmit apodization
    % use of SUM: summation over the number of delay series (e.g. MLT)
    DELAPOD = sum(exp(1i*kw*c*delaysTX),2).*APOD;
    
    %-- Summing the radiation patterns generating by all the elements
    RPk = RPmono*DELAPOD;
    %- include spectrum responses:
    RPk = pulseSPECT(k)*probeSPECT(k)*RPk;
    RPk(isOUT) = 0;
       
    %-- Output
    if isMKMOVIE % for MKMOVIE only (spectrum of the pressure field)
        SPECT(k,:) = RPk(1:nx);
        if ~isempty(options.RC)
            SPECT(k,:) = SPECT(k,:) +...
                (RPk(nx+1:end).*options.RC(:)).'*EXP_RC;                
        end
    elseif isSIMUS % Receive: for SIMUS only (spectra of the RF signals)
        if siz0(2) == 1
           SPECT(k,:) = probeSPECT(k) *... % the array bandwidth is considered
            (RPk.*options.RC(:)).'*RPmono ... % pressure received by the elements
            ; % *f(k)^2/fc^2; % Rayleigh scattering (OPTIONAL)
        else
            RPk = reshape( RPk.*options.RC(:), siz0 );
            RPmono = reshape( RPmono, [siz0 size( RPmono,2)] );
            SPECT(k,:,:) = probeSPECT(k) *... % the array bandwidth is considered
                squeeze( sum( RPk.*RPmono,1) ).' ... % pressure received by the elements
                ; % *f(k)^2/fc^2; % Rayleigh scattering (OPTIONAL)
        end
        if any(param.RXdelay) % reception delays, if any
            SPECT(k,:,:) = SPECT(k,:,:).*exp(1i*kw*c*param.RXdelay);
        end
    else % using PFIELD alone
        RP = RP + abs(RPk).^2; % acoustic intensity
        if nargout==3
            SPECT(k,:) = RPk;
        end
    end
    
    %- update the wait bar
    if options.WaitBar && rem(k,10)==0
        tstep = toc(tstart);
        trem = tstep*(nSampling/k-1);
        waitbar(k/nSampling,hwb,...
            {['Elapsed: ',int2str(floor(tstep/60)) ,' min ',...
            int2str(floor(rem(tstep,60))),' s'],...
            ['Remaining: ',int2str(floor(trem/60)) ,' min ',...
            int2str(floor(rem(trem,60))),' s']})
    end
    
end

%------------------------------------%
% end of SUMMATION OVER THE SPECTRUM %
%------------------------------------%



% Close the wait bar
if options.WaitBar, close(hwb), end

% Correcting factor (including integration step, df)
if isinf(NoW)
    CorFac = 1;
else
    CorFac = df;
end
if ~ElevationFocusing % 2-D space
    CorFac = CorFac*ElementWidth;
end
if exist('SPECT','var'), SPECT = SPECT*CorFac; end
RP = RP*CorFac;

% RMS acoustic pressure (if we are in PFIELD only)
if ~(isSIMUS || isMKMOVIE)
    RP = reshape(sqrt(RP),siz0);
    if nargout>2
        SPECT = reshape(SPECT,[siz0 nSampling]);
    end
end

end



function [P,param] = RunTheExample

% 2.7 MHz phased-array
param = getparam('P4-2v');

% Location of the focus
figure('Name','MUST - Matlab Ultrasound Toolbox - Damien Garcia',...
    'NumberTitle','off');
plot((param.Nelements-1)*param.pitch/2*[-1 1],[0 0],'k','Linewidth',10)
axis equal ij
axis([-4e-2 4e-2 0 10e-2])
xlabel('x-position (m)')
ylabel('z-position (m)')
hp = helpdlg({'We have a 2.7 MHz phased array.',...
    'Choose the focus point in the figure.'},'Focus point selection');
waitfor(hp)
[x0,z0] = ginput(1);
hold on
plot(x0,z0,'k*')
hold off

% TX time delays
dels = txdelay(x0,z0,param);

% Grid
x = linspace(-4e-2,4e-2,200);
z = linspace(0,10e-2,200);
[x,z] = meshgrid(x,z);

% RMS pressure field
opt.WaitBar = true;
P = pfield(x,zeros(size(x)),z,dels,param,opt);

% Final figure
imagesc(x(1,:)*1e2,z(:,1)*1e2,20*log10(P/max(abs(P(:)))))
caxis([-20 0])
c = colorbar;
c.YTickLabel{end} = '0 dB';
title('RMS pressure field')
colormap hot
hold on
plot((param.Nelements-1)*param.pitch/2*[-1 1]*1e2,param.pitch*[1 1]*1e2,...
    'g','Linewidth',6)
axis equal ij
axis([-4 4 0 10])
plot(x0*1e2,z0*1e2,'k*')
xlabel('[cm]')
hold off

end


function f = fresnelint(x)

% FRESNELINT Fresnel integral.
%
% J = FRESNELINT(X) returns the Fresnel integral J = C + 1i*S.
%
% We use the approximation introduced by Mielenz in
%       Klaus D. Mielenz, Computation of Fresnel Integrals. II
%       J. Res. Natl. Inst. Stand. Technol. 105, 589 (2000), pp 589-590
%

siz0 = size(x);
x = x(:);

issmall = abs(x)<=1.6;
c = zeros(size(x));
s = zeros(size(x));

% When |x| < 1.6, a Taylor series is used (see Mielenz's paper)
if any(issmall)
    n = 0:10;
    cn = [1 cumprod(-pi^2*(4*n+1)./(4*(2*n+1).*(2*n+2).*(4*n+5)))];
    sn = [1 cumprod(-pi^2*(4*n+3)./(4*(2*n+2).*(2*n+3).*(4*n+7)))]*pi/6;
    n = [n 11];
    c(issmall) = sum(cn.*x(issmall).^(4*n+1),2);
    s(issmall) = sum(sn.*x(issmall).^(4*n+3),2);    
end

% When |x| > 1.6, we use the following:
if any(~issmall)
    n = 0:11;
    fn = [0.318309844, 9.34626e-8, -0.09676631, 0.000606222, ...
        0.325539361, 0.325206461, -7.450551455, 32.20380908, ...
        -78.8035274, 118.5343352, -102.4339798, 39.06207702];
    gn = [0, 0.101321519, -4.07292e-5, -0.152068115, -0.046292605, ...
        1.622793598, -5.199186089, 7.477942354, -0.695291507, ...
        -15.10996796, 22.28401942, -10.89968491];
    fx = sum(fn.*x(~issmall).^(-2*n-1),2);
    gx = sum(gn.*x(~issmall).^(-2*n-1),2);    
    c(~issmall) = 0.5*sign(x(~issmall)) + ...
        fx.*sin(pi/2*x(~issmall).^2) - gx.*cos(pi/2*x(~issmall).^2);
    s(~issmall) = 0.5*sign(x(~issmall)) - ...
        fx.*cos(pi/2*x(~issmall).^2) - gx.*sin(pi/2*x(~issmall).^2);
end

f = reshape(c,siz0) + 1i*reshape(s,siz0);

end


function [A,B] = MGBMcoeff(n)

% Coefficients for the Multi-Gaussian beam models

switch n
    case 1
        A = 1.186;
        B = 3.92;
    case 2
        A = 0.452 - 0.875i;
        B = 5.021 - 6.61i;
        A = [A; conj(A)];
        B = [B; conj(B)];
    case 3
        A = -0.538 - 0.104i;
        B = 4.81 - 15.506i;
        A = [2.141 ; A ; conj(A)];
        B = [8.03 ; B ; conj(B)];
    case 4
        A = [0.187 + 0.275i ; 0.288 - 1.954i];
        B = [4.558 - 25.59i ; 8.598 - 7.924i];
        A = [A; conj(A)];
        B = [B; conj(B)];
    case 5
        A = [0.057 - 0.292i ; -1.79 + 0.385i];
        B = [5.078 - 34.993i ; 9.982 - 16.61i];
        A = [4.509 ; A ; conj(A)];
        B = [12.404 ; B ; conj(B)];
    case 10
        % Wen JJ, Breazeale MA. A diffraction beam field expressed as the
        % superposition of Gaussian beams. The Journal of the Acoustical
        % Society of America. 1988 May;83(5):1752-6.
        A = [11.4280 + 0.9517i ; 0.0600 - 0.0801i ; -4.2743 - 8.5562i ;...
            1.6576 + 2.7015i ; -5.0418 + 3.2488i ; 1.1227 - 0.6885i ;...
            -1.0106 - 0.2696i ; -2.5974 + 3.2202i ; -0.1484 - 0.3119i ;...
            -0.2085 - 0.2385i];
        B = 4*[4.0697 + 0.2273i ; 1.1531 - 20.9330i ; 4.4608 + 5.1268i ;...
            4.3521 + 14.9970i ; 4.5443 + 10.0030i ; 3.8478 + 20.0780i ;...
            2.5280 - 10.3100i ; 3.3197 - 4.8008i ; 1.9002 - 15.8200i ;...
            2.6340 + 25.0090i];
    case 10.5
        % Liu W, Ji P, Yang J. Development of a simple and accurate
        % approximation method for the Gaussian beam expansion technique.
        % Journal of the Acoustical Society of America. 2008
        % May;123(5):3516.
        A = [-0.0366 + 0.0698i ; -0.2880 - 0.1072i ; 0.0463 - 0.8593i ;...
            2.4278 - 0.4273i ; -1.6515 + 6.9321i];
        A = [A; conj(A(end:-1:1))];
        B = [0.9568 + 22.0499i ; 1.8966 + 17.3281i ; 2.5687 + 12.2845i ;...
            3.1522 + 7.1375i ; 3.7397 + 2.2497i];
        B = 4*[B; conj(B(end:-1:1))];
    case 15
        % Liu W, Ji P, Yang J. Development of a simple and accurate
        % approximation method for the Gaussian beam expansion technique.
        % Journal of the Acoustical Society of America. 2008
        % May;123(5):3516.
        A = [-0.0647 - 0.0042i ; 0.0334 - 0.2398i ; 0.5113 - 0.0972i ;...
            0.5858 + 0.7912i ; -0.6908 + 1.5627i ; -3.0363 + 0.5081i ;...
            -3.6501 - 6.3857i ; 13.6222];
        A = [A; conj(A(end-1:-1:1))];
        B = [1.2100 + 35.6867i ; 2.3108 + 31.3481i ; 2.8161 + 26.2901i ;...
            3.2223 + 21.1344i ; 3.4860 + 15.8696i ; 3.6537 + 10.4523i ;...
            4.0206 + 5.0002i ; 4.3552];
        B = 4*[B; conj(B(end-1:-1:1))];
    case 25
        % Liu W, Ji P, Yang J. Development of a simple and accurate
        % approximation method for the Gaussian beam expansion technique.
        % Journal of the Acoustical Society of America. 2008
        % May;123(5):3516.        
        A = [0.0250 + 0.0903i ; -0.1722 + 0.0417i ; -0.2335 - 0.2073i ;...
            0.0169 - 0.4669i ; 0.5003 - 0.4462i ; 0.9303 + 0.0207i ;...
            0.9359 + 0.8547i ; 0.2283 + 1.7047i ; -1.2466 + 2.0209i ;...
            -3.1897 + 1.1364i ; -5.3844 - 1.7978i ; -3.2923 - 12.5813i ;...
            22.7593];
        A = [A; conj(A(end-1:-1:1))];
        B = [2.7353 + 64.7628i ; 2.4159 + 60.8489i ; 3.0518 + 55.6123i ;...
            3.3108 + 50.1677i ; 3.5307 + 44.6995i ; 3.7217 + 39.1628i ;...
            3.8727 + 33.5837i ; 3.9930 + 27.9784i ; 4.0952 + 22.3454i ;...
            4.1578 + 16.6841i ; 4.2756 + 10.9041i ; 4.6409 + 5.2178i ;...
            4.9426];
        B = 4*[B; conj(B(end-1:-1:1))];
end

end


function structArray = IgnoreCaseInFieldNames(structArray)

switch inputname(1)
    case 'param'
        fieldLIST = {'attenuation','baffle','bandwidth','c','fc',...
            'fnumber','focus','fs','height','kerf','movie','Nelements',...
            'passive','pitch','radius','RXangle','RXdelay',...
            'TXapodization','TXfreqsweep','TXnow','t0','width'};
    case 'options'
        if isstruct(structArray)
            fieldLIST = {'dBThresh','ElementSplitting',...
                'FullFrequencyDirectivity','FrequencyStep','ParPool',...
                'WaitBar'};
        else
            return
        end
end

OldFieldNames = fieldnames(structArray);
tmp = lower(OldFieldNames);
assert(length(tmp)==length(unique(tmp)),...
    ['The structure ' upper(inputname(1)),...
    ' contains duplicate field names (when ignoring case).'])

[idx,loc] = ismember(lower(fieldLIST),tmp);
idx = find(idx); loc = loc(idx);
for k = 1:length(idx)
    tmp = eval(['structArray.' OldFieldNames{loc(k)}]); %#ok
    structArray = rmfield(structArray,OldFieldNames{loc(k)});
    eval(['structArray.' fieldLIST{idx(k)} ' = tmp;'])
end

end


function AdMessage
msgboxStruct.Interpreter = 'tex';
msgboxStruct.WindowStyle = 'modal';
try
[icondata,iconcmap] =...
    imread('http://www.biomecardio.com/images/MUSTicon.jpg');
catch
    icondata = ind2rgb(round(255*rescale(peaks(100))+1),...
        [1-hot(128); hot(128)]);
    iconcmap = [];
end
uiwait(msgbox(...
    {'\fontsize{9}MUST is an open-source toolbox for research purposes.',...
    'Please cite the corresponding articles.',...
    '',...
    'MUST \copyright 2020 Damien Garcia',...
    'LGPL-3.0-or-later','',...
    '\bf{garcia.damien@gmail.com}',...
    '\color[rgb]{0,0.5,0.5}\bullet \bf{www.biomecardio.com} \bullet'},...
    'MUST','custom',icondata,iconcmap,msgboxStruct));
end
