function [RF,param,RFspectrum] = simus(varargin)

%SIMUS   Simulation of ultrasound RF signals for a linear or convex array
%   RF = SIMUS(X,Y,Z,RC,DELAYS,PARAM) simulates ultrasound RF radio-
%   frequency signals generated by an ultrasound uniform LINEAR or CONVEX
%   array insonifying a medium of scatterers.
%   The scatterers are characterized by their coordinates (X,Y,Z) and
%   reflection coefficients RC.
%
%   X, Y, Z and RC must be of same size. The elements of the ULA are
%   excited at different time delays, given by the vector DELAYS. The
%   transmission and reception characteristics must be given in the
%   structure PARAM (see below for details).
%
%   RF = SIMUS(X,Z,RC,DELAYS,PARAM) or RF = SIMUS(X,[],Z,RC,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 "simus" in the command window for an example ---< 
%
%   The RF output matrix contains Number_of_Elements columns. Each column
%   therefore represents an RF signal. The number of rows depends on the
%   depth (estimated from max(Z)) and the sampling frequency PARAM.fs (see
%   below). By default, the sampling frequency is four times the center
%   frequency.
%
%   Units: X,Y,Z must be in m; DELAYS must be in s; RC has no unit.
%
%   DELAYS can also be a matrix. This alternative can be used to simulate
%   MLT (multi-line transmit) sequences. In this case, each ROW represents
%   a series of delays. For example, to create a 4-MLT sequence with a
%   64-element phased array, DELAYS matrix must have 4 rows and 64 columns
%   (size = [4 64]).
%
%   SIMUS uses PFIELD during transmission and reception. The parameters
%   that must be included in the structure PARAM are similar as those in
%   PFIELD. Additional parameters are also required (see below).
%
%   ---
%   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 (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" in "help pfield"
%
%       *** 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.
%
%       *** RECEIVE PARAMETERS *** (not in PFIELD)
%   14) PARAM.fs: sampling frequency (in Hz, default = 4*param.fc)
%   15) PARAM.RXdelay: reception law delays (in s, default = 0)
%
%   Other syntaxes:
%   --------------
%   i}  [RF,PARAM] = SIMUS(...) updates the fields of the PARAM structure.
%   ii} [...] = SIMUS without any input argument provides an interactive
%        example designed to produce RF signals from a focused ultrasound
%        beam using a 2.7 MHz phased-array transducer (without elevation
%        focusing).
%
%   PARALLEL COMPUTING:
%   ------------------
%   SIMUS calls the function PFIELD. If you have the Parallel Computing
%   Toolbox, SIMUS can execute several PFIELDs in parallel. If this option
%   is activated, a parallel pool is created on the default cluster. All
%   workers in the pool are used. The X and Z are splitted into NW chunks,
%   NW being the number of workers. To execute parallel computing, use: 
%       [...] = SIMUS(...,OPTIONS),
%   with OPTIONS.ParPool = true (default = false).
%
%
%   OTHER OPTIONS:
%   -------------
%      %-- FREQUENCY STEP & FREQUENCY SAMPLES --%
%   1a) Only frequency components of the transmitted signal in the range
%       [0,2fc] with significant amplitude are considered. The default
%       relative amplitude is -100 dB. You can change this value by using
%       the following:
%           [...] = SIMUS(...,OPTIONS),
%       where OPTIONS.dBThresh is the threshold in dB (default = -100).
%   1b) The frequency step is determined automatically to avoid aliasing in
%       the time domain. This step can be adjusted with a scaling factor
%       OPTIONS.FrequencyStep (default = 1). It is not recommended to
%       modify this scaling factor in SIMUS.
%   ---
%      %-- FULL-FREQUENCY DIRECTIVITY --%   
%   2)  By default, the directivity of the elements depends only on the
%       center frequency. This makes the calculation faster. To make the
%       directivities fully frequency-dependent, use: 
%           [...] = SIMUS(...,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) SIMUS calls the function PFIELD. It 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 a Fraunhofer
%      (far-field) approximation. 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 SIMUS 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 SIMUS will be submitted by February 2021.
%
%
%   A simple EXAMPLE:
%   ----------------
%   %-- Generate RF signals using a phased-array transducer
%   % Phased-array @ 2.7 MHz:
%   param = getparam('P4-2v');
%   % TX time delays:
%   x0 = 0; z0 = 3e-2; % focus location 
%   dels = txdelay(x0,z0,param);
%   % Six scatterers:
%   x = zeros(1,6); y = zeros(1,6);
%   z = linspace(1,10,6)*1e-2;
%   % Reflectivity coefficients:
%   RC = ones(1,6);
%   % RF signals:
%   param.fs = 20*param.fc; % sampling frequency
%   RF = simus(x,y,z,RC,dels,param);
%   % Plot the RF signals
%   plot(bsxfun(@plus,RF(:,1:7:64)/max(RF(:)),(1:10))',...
%      (0:size(RF,1)-1)/param.fs*1e6,'k')
%   set(gca,'XTick',1:10,'XTickLabel',int2str((1:7:64)'))
%   title('RF signals')
%   xlabel('Element number'), ylabel('time (\mus)')
%   xlim([0 11]), axis ij 
%   
%
%   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 PFIELD, TXDELAY, MKMOVIE, GETPARAM, GETPULSE.
%
%   -- Damien Garcia -- 2017/10, last update 2022/05/14
%   website: <a
%   href="matlab:web('https://www.biomecardio.com')">www.BiomeCardio.com</a>

if nargin==0
    if nargout>0
        [RF,param] = RunTheExample;
    else
        RunTheExample;
    end
    return
end

narginchk(5,7)
nargoutchk(0,3)

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

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


%------------------------%
% CHECK THE INPUT SYNTAX % 
%------------------------%


param = IgnoreCaseInFieldNames(param);
options = IgnoreCaseInFieldNames(options);
options.CallFun = 'simus';

%-- 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).')

%-- Parallel pool
if ~isfield(options,'ParPool')
    options.ParPool = false;
end

%-- Check if syntax errors may appear when using PFIELD
try
    opt = options;
    opt.ParPool = false;
    opt.WaitBar = false;
    [~,param] = pfield([],[],delaysTX,param,opt);
catch ME
    throw(ME)
end

%-- Sampling frequency (in Hz)
if ~isfield(param,'fs')
    param.fs = 4*param.fc; % default
end
assert(param.fs>=4*param.fc,'PARAM.fs must be >= 4*PARAM.fc.')

NumberOfElements = param.Nelements; % number of array elements

%-- Receive delays (in s)
if ~isfield(param,'RXdelay')
    param.RXdelay = zeros(1,NumberOfElements);
else
    assert(isvector(param.RXdelay) && isnumeric(param.RXdelay),...
        'PARAM.RXdelay must be a vector')
    assert(numel(param.RXdelay)==NumberOfElements,...
        'PARAM.RXdelay must be of length = (number of elements)')
    param.RXdelay = reshape(param.RXdelay,1,NumberOfElements);
end

%-- dB threshold (in dB: faster computation if lower value)
if ~isfield(options,'dBThresh')
    options.dBThresh = -100; % default is -100dB in SIMUS
end
assert(isscalar(options.dBThresh) && isnumeric(options.dBThresh) &&...
    options.dBThresh<0,'OPTIONS.dBThresh must be a negative scalar.')

%-- Frequency step (scaling factor)
% The frequency step is determined automatically. It is tuned to avoid
% aliasing in the temporal domain. The frequency step can be adjusted by
% using a scaling factor. 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.')
if options.FrequencyStep>1
    warning('MUST:FrequencyStep',...
            'OPTIONS.FrequencyStep is >1: aliasing may be present!')
end

RadiusOfCurvature = param.radius; % array radius-of-curvature (in m)


%-------------------------------%
% end of CHECK THE INPUT SYNTAX %
%-------------------------------%


%-- Centers of the tranducer elements (x- and z-coordinates)
if isinf(RadiusOfCurvature)
    % Linear array
    xe = ((0:NumberOfElements-1)-(NumberOfElements-1)/2)*param.pitch;
    ze = zeros(1,NumberOfElements);
else
    % 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

%-- Maximum distance
d2 = (x(:)-xe).^2+(z(:)-ze).^2;
maxD = max(sqrt(d2(:))); % maximum element-scatterer distance
clear d2
maxD = maxD + max(delaysTX(:))*param.c;

%-- Split x, y, z and RC if using a PARALLEL POOL of workers
assert(isscalar(options.ParPool) && islogical(options.ParPool),...
    'OPTIONS.ParPool must be a logical scalar (true or false).')
options.ParPool = options.ParPool &...
    license('test','Distrib_Computing_Toolbox');
if options.ParPool
    pool = gcp;
    NW = pool.NumWorkers; Nx = numel(x);
    dim1Dist = [ones(1,NW-1)*floor(Nx/NW) Nx-(NW-1)*floor(Nx/NW)];
    x = mat2cell(x(:),dim1Dist,1);
    if ElevationFocusing
        y = mat2cell(y(:),dim1Dist,1);
    end
    z = mat2cell(z(:),dim1Dist,1);
    RC = mat2cell(RC(:),dim1Dist,1);
end

%-- FREQUENCY SAMPLES
df = 1/2/(2*maxD/param.c); % to avoid aliasing in the time domain
df = df*options.FrequencyStep;
Nf = 2*ceil(param.fc/df)+1; % number of frequency samples


%-- Run PFIELD to calculate the RF spectra
RFspectrum = zeros(Nf,NumberOfElements, size(x,2) ); % will contain the RF spectra
options.FrequencyStep = df;
%- run PFIELD in a parallel pool (NW workers)
if options.ParPool
    options.WaitBar = false;
    spmd(NW)
        options.RC = RC{labindex};
        if ElevationFocusing
            [~,~,RFsp,idx] = pfield(x{labindex},y{labindex},z{labindex},...
                delaysTX,param,options);
        else
            [~,~,RFsp,idx] = pfield(x{labindex},z{labindex},...
                delaysTX,param,options);
        end
    end
    % sum the NW spectra (RFsp is a composite object)
    for k = 1:NW
        RFspectrum(idx{k},:) = RFspectrum(idx{k},:) + RFsp{k};
    end
else
    %- no parallel pool 
    options.RC =  RC;
    [~,~,RFsp,idx] = pfield(x,y,z,delaysTX,param,options);
    RFspectrum(idx,:,:) = RFsp;
end

%-- RF signals (in the time domain)
nf = ceil(param.fs/2/param.fc*(Nf-1));
RF = ifft(conj(RFspectrum),nf,'symmetric');
RF = RF(1:floor(nf/2),:,:)*param.fs/4/param.fc;

%-- Zeroing the very small values
RelThresh = 1e-5; % -100 dB
tmp = @(RelRF) 0.5*(1+tanh((RelRF-RelThresh)/(RelThresh/10)));
tmp = @(RelRF) round(tmp(RelRF)/(RelThresh/10))*(RelThresh/10);
RF = RF.*tmp(abs(RF)/max(abs(RF(:))));

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','TXdelay','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 [RF,param] = RunTheExample

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

% Location of the scatterers
Ns = 10; % number of scatterers
%
figure('Name','MUST - Matlab Ultrasound Toolbox - Damien Garcia',...
    'NumberTitle','off');
plot((param.Nelements-1)*param.pitch/2*[-1 1],[0 0],'g','Linewidth',5)
hold on
axis equal ij
axis([-4e-2 4e-2 0 10e-2])
xlabel('x-position (m)')
ylabel('z-position (m)')
hp = helpdlg({'We want a diverging wave with a 2.7 MHz phased array.',...
    ['Choose ' int2str(Ns) ' points (by clicking) in the figure.']},...
    'Scatterer selection');
waitfor(hp)
x = zeros(1,Ns); z = x;
for k = Ns:-1:1
    if k==Ns
        title(['Choose ' int2str(Ns) ' points'])
    elseif k==1
        title('A last one')
    else
        title(['Still ' int2str(k) ' points'])
    end
    [x(k),z(k)] = ginput(1);
    scatter(x(k),z(k),8,'r','filled')
end
hold off

% Add one null-RC point to force propagation till 10 cm
x = [x 0];
z = [z 10e-2];
RC = [ones(1,Ns) 0];

% TX time delays (diverging wave, angular width = 90 degrees):
dels = txdelay(param,0,pi/2);

% RF signals
param.fs = 10*param.fc; % sampling frequency
[RF,param] = simus(x,z,RC,dels,param);

% Plot the RF signals
plot((RF(:,1:3:64)/max(RF(:))*2+(1:3:64))',...
    (0:size(RF,1)-1)/param.fs*1e6,'k')
set(gca,'XTick',1:3:64,'XTickLabel',int2str((1:3:64)'))
title('RF signals')
xlabel('Element number'), ylabel('time (\mus)')
axis([-1 66 0 size(RF,1)/param.fs*1e6]), axis ij

ButtonName = questdlg('Do you want a movie with MKMOVIE?','MKOVIE',...
    'Yes','No','Yes');

switch ButtonName
    case 'Yes'
        % Image size (in cm)
        param.movie = [8 10];
        
        % Movie frames
        title('Generating the movie...')
        [F,info,param] = mkmovie(x,z,RC,dels,param);
        
        % Check the movie frames
        colormap([1-hot(128); hot(128)]);
        for k = 1:size(F,3)
            image(info.Xgrid*100,info.Zgrid*100,F(:,:,k))
            hold on
            plot((param.Nelements-1)*param.pitch/2*[-100 100],...
                [0 0],'g','Linewidth',5)
            scatter(x(1:10)*100,z(1:10)*100,8,'r','filled')
            hold off
            axis equal ij
            axis([-4 4 0 10])
            xlabel('x-position (cm)')
            ylabel('z-position (cm)')
            title([int2str(info.TimeStep*k*1e6) ' \mus'])
            drawnow
        end
        
    case 'No'
        % do nothing
    otherwise
        % do nothing
end

end


