function [M,T] = pp2bsp(pp)
% PP2BSP converts function in pp-form to genlib b-spline format
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

X = pp.breaks;
k = pp.order;
dim = pp.dim;
nX = numel(X);
dX = (X(end)-X(1))/(nX-1);
% Use triple knots to enforce only zero-order continuity
T = [X(1)+dX*(-k+1:0),reshape(repmat(X(2:end-1),k-1,1),1,[]),X(end)+dX*(0:k-1)];

% Fit based on k values in each segment
Z = (bsxfun(@times,pp.breaks(1:end-1),(1:k).'   )+...
     bsxfun(@times,pp.breaks(2:end  ),(k:-1:1).'))/(k+1);
Y = reshape(ppval(pp,Z(:)),prod(dim),numel(Z)).';
c = bspsum(T,eye(numel(T)-k),Z(:));
M = reshape(c\Y,[numel(T)-k,dim]);
