function [A,Psi]=RBF_QR_diffmat_2D(op,xe, varargin)
% Alternative calls:
% [A,Psi]=RBF_QR_diffmat(op,xe,xk,ep) % First time 
% [A,Psi]=RBF_QR_diffmat(op,xe,Psi)   % If Psi is already computed
%
% Computes a differentiation matrix (weights for RBF-FD stencils)
% using Gaussian RBFs at evaluation point(s) xe for RBFs placed at
% the nodes xk. 
%
%--- op (char)   : Alternatives '1', 'x', 'y', 'xx', 'xy', 'yy' ,'L',
%                   'Ln' where n is a number indicating the degree
%                   of the Laplacian.
%--- xe(1:M,1:2)  : The evaluation point(s) in Cartesian coordinates
%--- xk(1:N,1:2)  : The node points (no particular scaling assumed)
%--- ep (scalar)  : The (constant) shape parameter 
%--- Psi (struct) : Generated by a call to this function. Defines
%                   the RBF-QR basisfunctions Psi.
%
%--- Check if Psi exists. Otherwise, scale nodes and compute Psi.
if length(varargin)>=2
  xk = varargin{1};
  ep = varargin{2};
  %
  % If only one evaluation point, place that in zero.
  %
  if (size(xe,1)==1)
    cc = xe;
  else  
    cc = sum(xk,1)/size(xk,1)                    % Center of disc
  end  
  xk = [xk(:,1)-cc(1) xk(:,2)-cc(2)];
  xe = [xe(:,1)-cc(1) xe(:,2)-cc(2)];  
  r = sqrt(sum(xk.^2,2));
  re = sqrt(sum(xe.^2,2));
  rr = max(max(r),max(re));          % Radius of disk
    
  xk = [(1/rr)*r  atan2(xk(:,2),xk(:,1))]; % Polar coordinates
  xe = [(1/rr)*re atan2(xe(:,2),xe(:,1))]; % Polar coordinates
  ep = ep*rr;

  [Psi] = InitPsi_2D(ep,xk);
  Psi.rr = rr;
  Psi.cc = cc;

  %--- Also compute the interpolation matrix which is reused for all ops 
  Psi.A0 = RBF_QR_mat_2D(Psi,'1',Psi.xk);
  
elseif length(varargin)==1
  Psi = varargin{1};
  xe = [xe(:,1)-Psi.cc(1) xe(:,2)-Psi.cc(2)];
  re = sqrt(sum(xe.^2,2));
  xe = [(1/Psi.rr)*re atan2(xe(:,2),xe(:,1))]; % Polar coordinates  
end  

%--- Compute the differentiation matrix/ces
if (~iscellstr(op))
  %--- We are just computing one operator
  A = RBF_QR_mat_2D(Psi,op,xe)/Psi.A0;
  A = rescale_op(A,Psi.rr,op);
else
  numop = length(op);
  A = cell(1,numop);
  for i = 1:numop
    if (i==1) var=xe; else var=T; end
    op{i}
    Psi.ep
    Psi.xk
    xe
    [A{i},T] = RBF_QR_mat_2D(Psi,op{i},var);
    A{i} = A{i}/Psi.A0;
    A{i} = rescale_op(A{i},Psi.rr,op{i});
  end
end

function A = rescale_op(A,rr,op)
%--- Adjust the scaling according to derivative
[deg,diff,op]=RBF_QR_parse(op);
A = (1/rr).^deg*A;

