pd = 0.9;            % Probability of detection
pfa = 1e-6;          % Probability of false alarm
max_range = 5000;    % Maximum unambiguous range
tgt_rcs = 1;         % Required target radar cross section
int_pulsenum = 10;   % Number of pulses to integrate

load BasicMonostaticRadarExampleData;

fc = radiator.OperatingFrequency;         % Operating frequency (Hz)
v = radiator.PropagationSpeed;            % Wave propagation speed (m/s)
lambda = v/fc;                            % Wavelength (m)
fs = waveform.SampleRate;                 % Sampling frequency (Hz)
prf = waveform.PRF;                       % Pulse repetition frequency (Hz)

ura = phased.URA('Element',antenna,...
    'Size',[30 30],'ElementSpacing',[lambda/2, lambda/2]);
% Configure the antenna elements such that they only transmit forward
ura.Element.BackBaffled = true;

radiator.Sensor = ura;
collector.Sensor = ura;

% We need to set the WeightsInputPort property to true to enable it to
% accept transmit beamforming weights
radiator.WeightsInputPort = true;

% Calculate the array gain
arraygain = phased.ArrayGain('SensorArray',ura,'PropagationSpeed',v);
ag = arraygain(fc,[0;0]);

% Use the radar equation to calculate the peak power
snr_min = albersheim(pd, pfa, int_pulsenum);
peak_power = radareqpow(lambda,max_range,snr_min,waveform.PulseWidth,...
    'RCS',tgt_rcs,'Gain',transmitter.Gain + ag);

% Set the peak power of the transmitter
transmitter.PeakPower = peak_power;

initialAz = 90; endAz = -90;
volumnAz = initialAz - endAz;

% Calculate 3-dB beamwidth
theta = radtodeg(sqrt(4*pi/db2pow(ag)));

scanstep = -6;
scangrid = initialAz+scanstep/2:scanstep:endAz;
numscans = length(scangrid);
pulsenum = int_pulsenum*numscans;

% Calculate revisit time
revisitTime = pulsenum/prf;

tgtpos = [[2800; 3200; 0],[2500; 600; 0]];
tgtvel = [[0; 0; 0],[0; 0; 0]];
tgtmotion = phased.Platform('InitialPosition',tgtpos,'Velocity',tgtvel);

tgtrcs = [1.6 2.2];
target = phased.RadarTarget('MeanRCS',tgtrcs,'OperatingFrequency',fc);

% Calculate the range, angle, and speed of the targets
[tgtrng,tgtang] = rangeangle(tgtmotion.InitialPosition,...
    sensormotion.InitialPosition);

numtargets = length(target.MeanRCS);

% Create the steering vector for transmit beamforming
steeringvec = phased.SteeringVector('SensorArray',ura,...
    'PropagationSpeed',v);

% Create the receiving beamformer
beamformer = phased.PhaseShiftBeamformer('SensorArray',ura,...
    'OperatingFrequency',fc,'PropagationSpeed',v,...
    'DirectionSource','Input port');

% Define propagation channel for each target
channel = phased.FreeSpace(...
    'SampleRate',fs,...
    'TwoWayPropagation',true,...
    'OperatingFrequency',fc);

fast_time_grid = unigrid(0, 1/fs, 1/prf, '[)');

% Pre-allocate array for improved processing speed
rxpulses = zeros(numel(fast_time_grid),pulsenum);

for m = 1:pulsenum

    % Update sensor and target positions
    [sensorpos,sensorvel] = sensormotion(1/prf);
    [tgtpos,tgtvel] = tgtmotion(1/prf);

    % Calculate the target angles as seen by the sensor
    [tgtrng,tgtang] = rangeangle(tgtpos,sensorpos);

    % Calculate steering vector for current scan angle
    scanid = floor((m-1)/int_pulsenum) + 1;
    sv = steeringvec(fc,scangrid(scanid));
    w = conj(sv);

    % Form transmit beam for this scan angle and simulate propagation
    pulse = waveform();
    [txsig,txstatus] = transmitter(pulse);
    txsig = radiator(txsig,tgtang,w);
    txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel);

    % Reflect pulse off of targets
    tgtsig = target(txsig);

    % Beamform the target returns received at the sensor
    rxsig = collector(tgtsig,tgtang);
    rxsig = receiver(rxsig,~(txstatus>0));
    rxpulses(:,m) = beamformer(rxsig,[scangrid(scanid);0]);
end

% Matched filtering
matchingcoeff = getMatchedFilter(waveform);
matchedfilter = phased.MatchedFilter(...
    'Coefficients',matchingcoeff,...
    'GainOutputPort',true);
[mf_pulses, mfgain] = matchedfilter(rxpulses);
mf_pulses = reshape(mf_pulses,[],int_pulsenum,numscans);

matchingdelay = size(matchingcoeff,1)-1;
sz_mfpulses = size(mf_pulses);
mf_pulses = [mf_pulses(matchingdelay+1:end) zeros(1,matchingdelay)];
mf_pulses = reshape(mf_pulses,sz_mfpulses);

% Pulse integration
int_pulses = pulsint(mf_pulses,'noncoherent');
int_pulses = squeeze(int_pulses);

% Visualize
r = v*fast_time_grid/2;
X = r'*cosd(scangrid); Y = r'*sind(scangrid);
clf;
radar2dplot = pcolor(X,Y,pow2db(abs(int_pulses).^2));
axis equal tight
shading interp
axis off

rotate(radar2dplot, [ 0 0 1], 90 );