function L = measureLobes(angd,mag, units)
%measureLobes Measure antenna lobes and metrics from directivity pattern.
%   L = measureLobes(A,M) returns main, back and side lobe measurements
%   from the antenna directivity pattern described by angle vector A and
%   magnitude vector M, assuming polar data with magnitude in dB gain and
%   angle in degrees.
%
%   If omitted, A is taken to be uniformly distributed from 0 to 360
%   degrees, where A = (0:N-1)/N*360 and N = numel(M).
%
%   L is a structure with antenna lobe measurements and antenna
%   performance measurements in each of its fields.
%
%   ANTENNA LOBE MEASUREMENTS
%     L.main contains a structure describing the main lobe.
%       L.main.magnitude
%       L.main.angle
%       L.main.index
%       L.main.extent
%
%     L.back contains a structure describing the back lobe.
%       L.back.magnitude
%       L.back.angle
%       L.back.extent
%
%     L.side contains a structure describing the side lobes.
%       L.side.magnitude
%       L.side.angle
%       L.side.index
%       L.side.extent
%
%   ANTENNA PERFORMANCE METRICS
%     L.FB: Front-to-back ratio (dB)
%        Maximum magnitude of main lobe minus the interpolated magnitude
%        level 180-degrees from the maximum
%
%     L.FBIdx: Front-to-back ratio peak indices
%        1x2 vector of indices into A and M at which the front- and
%        back-lobe peaks were measured for Front-to-Back Ratio
%
%     L.SLL: Sidelobe level relative to main lobe (dB)
%        Maximum of main lobe minus highest side- or back-lobe level.
%
%     L.SLLIdx: Sidelobe level magnitude indices
%        1x2 vector of indices into A and M at which the front- and side-
%        lobe peaks were measured
%
%     L.FNBW: First-null beamwidth (degrees)
%        Angle between the first nulls appearing on each side of the main
%        lobe peak
%
%     L.FNBWIdx: First-null indices
%        1x2 vector of indices into A and M corresponding to the first
%        nulls on each side of the main lobe peak
%
%     L.HPBW: Half-power beamwidth (degrees)
%        Angle between magnitude points approximately half the power of the
%        maximum magnitude of the main lobe
%
%     L.HPBWIdx: Half-power beamwidth indices
%        1x2 vector of interpolated (fractional) indices at HPBW points
%
% See also findpeaks.

% Requires Signal Processing Toolbox due to use of findpeaks.

% Index vectors are in ascending order, to within angle-wraparound.

% Construct default return struct
%
% Default lobe info struct
L.mainLobe.index = [];
L.mainLobe.magnitude = [];
L.mainLobe.angle = [];
L.mainLobe.extent = []; % 1x2 index vector

% There's no index for back lobe value.
%
% We interpolate to get the angle and value at the back-reference angle
% (180-deg from main lobe peak).  It doesn't generally fall on a sample.
L.backLobe.magnitude = [];
L.backLobe.angle = [];
L.backLobe.extent = [];

% Only the highest sidelobe peak is recorded in index, mag and angle.  All
% side lobes are defined in Extent, however.
L.sideLobes.index = [];
L.sideLobes.magnitude = [];  % highest of the sidelobes
L.sideLobes.angle = [];
L.sideLobes.extent = [];  % all sidelobes

L.FB = [];
L.SLL = [];
L.HPBW = [];
L.FNBW = [];
L.FBIdx = [];
L.SLLIdx = [];
L.HPBWIdx = [];
L.HPBWAng = [];
L.FNBWIdx = [];
dbdown = 3;
interp_HPBW = 0; % Interpolation (linear) set to false
% ---

% Determine span of dataset
[angleSpan,fullCircle] = internal.polariCommon.findAngleSpan( ...
    angd*pi/180);

% Find the peaks:
%   findPolarPeaks(mag,circ,Npeaks,magLim,userOpts)
%   peaksIdx = internal.polariCommon.findPolarPeaks(mag,fullCircle,inf);
%   To enable wrapping of data circ is defined as true always
peaksIdx = internal.polariCommon.findPolarPeaks(mag,true,inf);

% Find the valleys:
valleysIdx = internal.polariCommon.findPolarValleys(mag,peaksIdx,fullCircle,units);

% MAIN LOBE
%
% At least one peak present
[L.mainLobe.magnitude,i] = max(mag(peaksIdx));
mainLobePeakIdx = peaksIdx(i);
peaksIdx(i) = []; % remove main peak from peak lists
L.mainLobe.index = mainLobePeakIdx;
L.mainLobe.angle = angd(mainLobePeakIdx);
Nv = numel(valleysIdx);
if Nv < 2
    % Cannot identify start/end of main lobe.
    %
    % *All* data is main lobe:
    L.mainLobe.extent = [1 numel(mag)];
    flag = ~any(isnan(mag))&& ~fullCircle;
    
    % Half-power Beamwidth (.HPBW)
    if flag 
        [L.HPBW, L.HPBWAng, L.HPBWIdx] = internal.calc_bw([NaN mag' NaN], 0,...
            [NaN angd' NaN], dbdown, interp_HPBW, 1, fullCircle);
        L.HPBWIdx(1) = find(angd == L.HPBWAng(1));
        L.HPBWIdx(2) = find(angd == L.HPBWAng(2));
    else
        [L.HPBW, L.HPBWAng, L.HPBWIdx] = internal.calc_bw(mag, 0,...
            angd, dbdown, interp_HPBW, 1, fullCircle);
    end
    return % EARLY EXIT
end

%  Finite extent of mainlobe:
% - closest valley index less than and greater than main lobe peak index
% - if we didn't find the index, it's because the main lobe fell before the
%   first valley identified at START of dataset (or fell after the last
%   valley identified at END of dataset).  The data is around the unit
%   circle, so that valley is found by wrapping to the end/start of the
%   data. 
%
valleysIdx = sort(valleysIdx);
i = find(valleysIdx < mainLobePeakIdx,1,'last');
if isempty(i)
    iLo = valleysIdx(end);
else
    iLo = valleysIdx(i);
end
i = find(valleysIdx > mainLobePeakIdx,1,'first');
if isempty(i)
    iHi = valleysIdx(1);
else
    iHi = valleysIdx(i);
end
% Arrange in CCW order
L.mainLobe.extent = [iLo iHi]; % ordered set

% SIDE LOBES
%
% All valley-to-valley intervals excluding the main lobe. We
% will identify and remove the back lobe from this list in a
% later step.

% Sidelobe extents are ordered pairs of adjacent valley
% indices including wraparound.  If Nv=4, extents is:
%     [1 2; 2 3; 3 4; 4 1],
% prior to removing main and back lobe extent rows.
%
% valleysIdx is already sorted.
% create Nx2 matrix of all valley index row-pairs:
idx = (1:Nv)';
sideLobeExtents = valleysIdx([idx circshift(idx,-1)]);

% Find row with main lobe indices
% - could subtract main lobe row from all rows and look for
%   [0 0]
% - simpler: find 1st index of extent pair matching 1st index
%  of main lobe extent
%bsxfun(@minus,sideLobeExtent,L.mainLobe.Extent)
idx = find(sideLobeExtents(:,1) - L.mainLobe.extent(1) == 0);
assert(~isempty(idx));
% Remove main lobe row, leaving side lobes + back lobe
sideLobeExtents(idx,:) = [];

% BACK LOBE
%
% Find lobe with extent that includes 180-deg away from main
% lobe.
%
% Also find its peak, which may not be exactly 180-deg away
% from main lobe peak.
%
% Start by computing the angle diametrically opposed to the
% peak of the main lobe:
mainAng = L.mainLobe.angle;
backRefAng = mainAng-180; % this angle may wrap

% Find lobe whose span includes backRefAng
%   -> that's the "back lobe"
%
% Does data span a full 360-deg?
% If not, check if backref angle is in the data span:
if fullCircle || ...
        internal.polariCommon.anglesWithinSpan(backRefAng,angleSpan)
    % backRefAng is within the data set -> necessary but
    % not sufficient to determine if we have a back lobe.
    %
    % No lobe may contain this angle, in which case there would
    % be no back lobe.
    
    % Need sideLobeAng for back-lobe determination:
    %
    % NOTE: indexing takes the shape of ang in general.
    % When sideLobeExtents is a matrix, the syntax happens to
    % return the shape of sideLobeExtents.  We want the shape
    % of sideLobeExtents always, which is Nx2.
    %
    tempLobeAng = angd(sideLobeExtents); % endpoints, degrees
    if isvector(tempLobeAng)
        tempLobeAng = tempLobeAng(:).'; % force into a row (1x2)
    end
    idx = internal.polariCommon.spanContainingAngle( ...
        tempLobeAng*pi/180, backRefAng*pi/180);
    if isempty(idx)
        % No back lobe (no lobe includes backrefang)
        backLobeExtent = [];
    else
        % Extract back lobe info, remove from side lobe list:
        assert(isscalar(idx))
        backLobeExtent = sideLobeExtents(idx,:);
        %backLobeExtAng = tempLobeAng(idx,:); % = angd(backLobeExtent)
        L.backLobe.extent = backLobeExtent;
        
        sideLobeExtents(idx,:) = []; % remove backlobe row
    end
else
    % Angle span of data does not include the angle
    % diametrically opposing the main lobe peak
    backLobeExtent = [];
end

% At this point, we have main lobe fully defined, as well as
% extents of back lobe and side lobes.
%
% Now we identify peaks of backlobe and each sidelobe.

% Find backlobe level, remove peaks in back lobe
%
if ~isempty(backLobeExtent)
    % Remove all peaks falling within back lobe
    %
    % backLobeExtAng: vector with two angles (a "span")
    % peaksIdx: vector of remaining peak indices
    %           (main lobe peak already removed)
    %{
                allPeakAng = angd(peaksIdx);
                sel = internal.polariCommon.anglesWithinSpan( ...
                    allPeakAng*pi/180, backLobeExtAng*pi/180);
                if any(sel)
                    % Only proceed if peaks found within back lobe extent
                    %idx = peaksIdx(sel);
                    % Remove backlobe peak from list of remaining peaks
                    peaksIdx(sel) = [];
                    
                    % Record peak value in backlobe
                    %
                    %L.backLobe.Index = idx;
                    %L.backLobe.Magnitude = mag(idx);
                    %L.backLobe.Angle = angd(idx);
                end
    %}
    
    % Record the value at backRefAng, whether it's a peak or
    % not.  That's the "back" level for front/back (F/B) ratio.
    %
    L.backLobe.angle = backRefAng;
    [L.backLobe.magnitude, L.backLobe.index] = ...
        internal.polariCommon.polarInterp(backRefAng,angd,mag);
end

% Sidelobes
%
% Record extent of all sidelobes, not including back lobe
L.sideLobes.extent = sideLobeExtents;
%
% Record the highest peak, not including main peak.
% This includes side lobes AND any back-lobe peaks.
%
[val,idx] = max(mag(peaksIdx));
L.sideLobes.index = peaksIdx(idx);
L.sideLobes.angle = angd(peaksIdx(idx));
L.sideLobes.magnitude = val;

% Compute Front/Back Ratio, Side Lobe Level, Beamwidth
%
m1 = L.mainLobe.magnitude;
if ~isempty(m1)
    % F/B Ratio (.FB)
    %
    % NOTE: FB is undefined if backlobe not identified
    %
    m2 = L.backLobe.magnitude;
    if ~isempty(m2)
        L.FB = m1 - m2;
        L.FBIdx = [L.mainLobe.index L.backLobe.index];
    end
    
    % Sidelobe Level (.SLL)
    %
    % NOTE: SLL is undefined if sidelobes not identified
    m2 = L.sideLobes.magnitude;
    if ~isempty(m2)
        L.SLL = m1 - m2;
        
        % xxx these indices may be in non-minimum span order
        %     (still correct, but possibly non-minimum span)
        L.SLLIdx = [L.mainLobe.index L.sideLobes.index];
    end
    
    % First-null Beamwidth (.FNBW)
    ext = L.mainLobe.extent;
    L.FNBW = internal.polariCommon.angleDiff(angd(ext)*pi/180)*180/pi;
    L.FNBWIdx = ext; % just a copy of mainlobe extent indices
    
    % Half-power Beamwidth (.HPBW)    
    [L.HPBW, L.HPBWAng, L.HPBWIdx] = internal.calc_bw(mag, 0, angd,...
            dbdown, interp_HPBW, 1, fullCircle);
end

end

%{
function [hpAng,closestIdx] = findHalfPowerAngles(L,angd,mag,fullCircle,dbdown)
% Compute half-power beamwidth, assuming mag is in dB.
%
% Assumes lobe info passed in, since we're in the midst of computing it.
%
% Returns closest indices corresponding to 3dB points, for visualization.

hpAng = [];
closestIdx = [];

% Get mag data over main lobe
ext = L.mainLobe.extent;
e1Idx = ext(1);
e2Idx = ext(2);

% Find 3dB-down values
magPk = L.mainLobe.magnitude;
if nargin < 5
    magdBdown = magPk-3;
else
    magdBdown = magPk-dbdown;
end
pkIdx = L.mainLobe.index;

% Find 1st 3dB point
%  - in range e1Idx to pkIdx
%  - given: e1Idx <= pkIdx
%  - m will be predominantly increasing from 1:end
if e1Idx > pkIdx
    if fullCircle
        m = mag([e1Idx:end 1:pkIdx]);
        a = angd([e1Idx:end 1:pkIdx]);
    else
        m = mag(e1Idx:end);
        a = angd(e1Idx:end);
    end
else
    m = mag(e1Idx:pkIdx);
    a = angd(e1Idx:pkIdx);
end
% Magnitude values may not be strictly monotonic increasing.
% interp1 will throw an error
% a1 = interp1(m,a,mag3dB);
[a1,i1] = interpLargelyIncrX(m,a,magdBdown);
if isempty(a1) % couldn't find point
    return
end

% Fractional interp index needs e1Idx offset and could wrap.
% Fix it up to be valid for vis purposes:
Nd = numel(mag);
i1 = i1+e1Idx-1; % compensate for index offset
if i1 > Nd
    i1 = i1-Nd+1;
end

% Find 2nd 3dB point
%  - in range pkIdx to e2Idx
%  - given: pkIdx >= e2Idx
%  - m will be predominantly increasing from 1:end
if pkIdx > e2Idx
    if fullCircle
        m = flip(mag([pkIdx:end 1:e2Idx]));
        a = flip(angd([pkIdx:end 1:e2Idx]));
    else
        m = flip(mag(pkIdx:end));
        a = flip(angd(pkIdx:end));
    end
else
    m = flip(mag(pkIdx:e2Idx));
    a = flip(angd(pkIdx:e2Idx));
end
[a2,i2] = interpLargelyIncrX(m,a,magdBdown);
if isempty(a2) % couldn't find point
    return
end

% Fractional interp index needs pkIdx offset and could wrap.
% Fix it up to be valid for vis purposes:
Nm = numel(m);          % compensate for "flip"
i2 = (Nm-i2+1)+pkIdx-1; % compensate for index offset
if i2 > Nd
    i2 = i2-Nd+1;
end

hpAng = [a1 a2];
closestIdx = [i1 i2];

end

function [av,idx] = interpLargelyIncrX(m,a,mv)
% Find interpolation point given x is "largely increasing"
%  - m may not be monotonic strictly increasing
%  - there may be more than one point that yields mv
%  - returns interp point closest to end of m (returns 'last')
%  - optionally returns fractional interp index
%  - a is the angle vector in degrees

% Find a pair of adjacent indices that bracket mv
%
% Find all points m < mv,
% removing i1==Nm since i2 would be out-of-bounds
i1 = find(m <= mv);
i1(i1 == numel(m)) = [];
if isempty(i1)
    av = [];
    idx = [];
    return
end
i1 = i1(end);  % Keep index closest to end of x

% Find value pair and interpolated value:
m1 = m(i1);
m2 = m(i1+1);
frac = (mv-m1)/(m2-m1);
idx = i1+frac;
av = a(i1) + frac * ...
    internal.polariCommon.angleDiffRel(a(i1:i1+1)*pi/180)*180/pi;

%{
% We *could* hit a pair of angles at exactly the wraparound point in
% degrees around the circle.  In that case, difference >= 180.
%
% It could be normal decrementing angle.  In that case, diff < 180.
%
if a_i1 > a_i2
    isWrap = internal.polariCommon.angleDiff >= 180;
    if isWrap
        % Wraparound at 0/360
        Nr = 1+floor((a_i1-a_i2)/360);
        a_i2 = a_i2+360*Nr; % a_i2 just past a_i1 modulo 360
    else
        % Decrementing angles
        t = a_i2; a_i2 = a_i1; a_i1 = t;
        frac = 1-frac;
    end
end

av = (1-frac)*a_i1 + frac*a_i2;
idx = i1+frac;
%}

end
%}
% [EOF]
