function [A,Psi]=RBF_QR_diffmat_1D(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 (chars)   : Alternatives '0', 'x', 'xx' or {'0','x',...}
%--- xe(1:Ne)  : The evaluation point(s) in Cartesian coordinates
%--- xk(1:N)  : 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};
  
  %--- Scale the nodes so that -1 and 1 are the extreme points
  x1 = min(min(xk),min(xe));
  x2 = max(max(xk),max(xe));
  cc = 0.5*(x1+x2);         % Center of points
  rr = 0.5*(x2-x1);         % radius (half length)
  
  xk = xk-cc;
  xe = xe-cc;

  xk = (1/rr)*xk;
  xe = (1/rr)*xe;
  ep = ep*rr;
  
  [Psi] = InitPsi_1D(ep,xk);
  Psi.rr = rr;
  Psi.cc = cc;

  %--- Also compute the interpolation matrix which is reused for all ops 
  Psi.A0 = RBF_QR_mat_1D(Psi,'0',Psi.xk);
  
elseif length(varargin)==1
  Psi = varargin{1};
  xe = xe-Psi.cc;
  xe = (1/Psi.rr)*xe;
end  

%--- Compute differentiation matrices for one or more operators 
if (~iscellstr(op))
  A=eval_diffmat(Psi,op,xe);
else
  for i=1:length(op)
    A{i}=eval_diffmat(Psi,op{i},xe);
  end
end


function A=eval_diffmat(Psi,op,xe)

%--- Compute the differentiation matrix
 A = RBF_QR_mat_1D(Psi,op,xe)/Psi.A0;

%--- Adjust the scaling according to derivative
if (strcmp(op,'x'))
  A=(1/Psi.rr)*A; 
elseif (strcmp(op,'xx'))
  A=(1/Psi.rr).^2*A;
end


