function meqplotarrows(rx,zx,Fx,Brx,Bzx,varargin)
%MEQPLOTARROWS: Flux surface plot with arrows indicating B field direction
% meqplotarrows(rx,zx,Fx,Brx,Bzx,varargin)
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

p=inputParser;
p.addParameter('parent',gca,@(x) ishandle(x));
p.addParameter('LineStyle','-');
p.addParameter('LineWidth',1);
p.addParameter('Color','k');
p.addParameter('ArrowLength',0.015); % arrow length [m]
p.addParameter('dec',10); % decimation of grid points (more=fewer arrows)
p.addParameter('plotContour',true); % Plot contour lines or not
p.addParameter('FxLevels',21); % contour values
p.KeepUnmatched = true;
parse(p,varargin{:}); P = p.Results;

hax = P.parent;
F = P.FxLevels;

if P.plotContour
  C = contour(rx,zx,Fx,F,'LineStyle',P.LineStyle,'Color',P.Color,'LineWidth',P.LineWidth,'parent',hax);
else
  C = contourc(rx,zx,Fx,F);
end
axis(hax,'equal'); axis(hax,'tight')

% interpolant
BBrx = griddedInterpolant({zx,rx},Brx);
BBzx = griddedInterpolant({zx,rx},Bzx);

len = P.ArrowLength;
ic=1;
while ic<size(C,2)
  nc = C(2,ic);
  rc = C(1,ic+1:P.dec:ic+nc);
  zc = C(2,ic+1:P.dec:ic+nc);
  ic = ic+nc+1;
 
  % interpolate
  Brc = BBrx(zc',rc')';
  Bzc = BBzx(zc',rc')';
  
  alpha = atan2(Brc,Bzc);
  ax	= rc - len*sin(pi/4+alpha);
  ay	= zc - len*cos(pi/4+alpha);
  bx	= rc + len*sin(pi/4-alpha);
  by	= zc - len*cos(pi/4-alpha);
  
  for ii=1:numel(rc)
    xx=[ax;rc;bx]; yy=[ay;zc;by];
    line(xx(:,ii),yy(:,ii),'Tag','arrows',...
      'LineWidth',P.LineWidth,'Color','k','parent',hax);
  end
end

end
