%LOCQ  Location of q surfaces
% [AR] = LOCRMEX(AQ,QQ,QR,NAR,NF) returns the location of the flux surfaces where
% QR=QQ(AQ==AR)
%
% AQ: x data (can be a matrix of size [nO,nQ] when QQ has nQ elements)
% QQ: y data
% QR: Requested y values (will be of size [nO,nR,naR] )
% NAR: (optional) Number of maximum R values to be found (per QR) for nonmonotonic profiles.
%          (default=1: outermost value found)
% NF: (optional) AR = NF signals the q value was not found
%
% See also LOCRMEX, LOCSMEX
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function aR = locQ(aQ,iqQ,iqR,naR,NF)

if nargin < 5, NF = NaN; end % default
if nargin < 4, naR = 1; end % default
if naR<1, error('invalid naR<1'); end
nR = numel(iqR);
nQ = numel(iqQ);
nO = size(aQ,1);

assert(size(aQ,2)==nQ,...
  'aQ must have numel(iqQ) columns')
aR = repmat(NF,nO,nR,naR); % default value for non-found points

for kR = 1:nR
  iq = iqR(kR);
  nrq = 1; % number of points found for this iq
  for kQ = nQ:-1:2
    if nrq>naR; break; end
    if (iq - iqQ(kQ)) * (iq - iqQ(kQ-1)) < 0 % interval contains point to be interpolated
      % linear interpolation
      aR(:,kR,nrq) = aQ(:,kQ) + (iqR(kR)-iqQ(kQ))*(aQ(:,kQ)-aQ(:,kQ-1))/(iqQ(kQ)-iqQ(kQ-1));
      nrq = nrq+1;
    elseif (iq==iqQ(kQ)) % this point
      aR(:,kR,nrq) = aQ(:,kQ);
      nrq = nrq+1;
    elseif (kQ==2) && (iq==iqQ(1)) % last point
      aR(:,kR,nrq) = aQ(:,1);    
      nrq = nrq+1;
    end
  end
end
end
