function qpM = qintc(inp,dr,dz)
% Auxiliary function for qint.m
% Pre-computes matrix to get coefficients for quadratic interpolation on N
% points. See qintmex.m for details.
%
% INP: Structure of parameters containing:
%   n: Number of points for interpolation (minimum 6, maximum 9)
%      if INP is empty, uses N=9 (default);
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% Returns ifp interpolation function parameter structure 
dr2=dr*dr; dz2=dz*dz; drdz=dr*dz;

if isempty(inp)
  n = 9;
else
  n = inp.n;
end

A = ...
 [1  0    0    0     0    0
  1  dr   0    0    dr2   0
  1 -dr   0    0    dr2   0
  1  0   dz    0     0   dz2
  1  0  -dz    0     0   dz2
  1  dr  dz  2*drdz dr2  dz2
  1 -dr -dz  2*drdz dr2  dz2
  1 -dr  dz -2*drdz dr2  dz2
  1  dr -dz -2*drdz dr2  dz2];
 
% Return pseudoinverse (transposed)
qpM = pinv(A(1:n,:)).';
