function [pp,ppd] = pprace(rk,zk,lk);

% PPRACE   Piecewise polynomial descrption of a racetrack
%   PPRACE(RK,ZK,LK) returns the piecewise polynomial PP describing the
%   racetrack given by its knots [RK,ZK]. PP is a complex valued function giving
%   R+iZ. LK specifies the length of the curved segment between knot k and k+1,
%   NaN for a straight segment. Optionaly returns the polynomial derivative.
%   Example for a circle: PPRACE([1 0 -1 0],[0 1 0 -1],pi/2*[1 1 1 1])
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if nargin == 0
 rk = [1 0 -1 0]; zk = [0 1 0 -1]; lk = pi/2*[1 1 1 1];
end

rk = rk(:); zk = zk(:); lk = lk(:); nk = length(rk);
rzk = rk + i*zk; drzk = diff([rzk;rzk(1)]);
kss = isnan(lk);

A = zeros(4*nk); B = zeros(4*nk,1);
for kk = 1:nk                                       
 k = (-3:0)+4*kk;                                  
 if kss(kk) 
  lk(kk) = abs(drzk(kk));
  % polynomial: a*s^3 + b*s^2 + c*s + d
  % a(k) = 0; b(k) = 0; c(k)*ds(k) = dx(k); d(k) = x(k);                               
  A(k,k) = diag([1 ; 1 ;   lk(kk) ; 1]);      
  B(k)   =      [0 ; 0 ; drzk(kk) ; rzk(kk)]; 
 else
  l = mod((-8:3)+4*kk,4*nk)+1;
  kk1 = mod(kk-2,nk)+1;
  % continuity at start:       d(k) = x(k)
  % continuity at end:         a(k)*l(k)^3 + b(k)*l(k)^2 + c(k)*l(k) + d(k) = x(k+1)
  % slope continuity at start: c(k) = 3*a(k-1)*l(k-1)^2 + 2*b(k-1)*l(k-1) + c(k-1)
  % slope continuity at end:   3*a(k)*l(k)^2 + b(k)*l(k) + c(k) = c(k+1)
  A(k,l) = A(k,l) + ...
  [0            0          0    0 0          0        0      1    0 0  0 0;
   0            0          0    0 lk(kk)^3   lk(kk)^2 lk(kk) 1    0 0  0 0;
  -3*lk(kk1)^2 -2*lk(kk1) -1    0 0          0        1      0    0 0  0 0;
   0            0          0    0 3*lk(kk)^2 2*lk(kk) 1      0    0 0 -1 0];
  B(k) = B(k) + [rzk(kk) ; rzk(mod(kk,nk)+1) ; 0 ; 0];
  if ~kss(kk1) % two contiguous curved segments
   % curvature continuity at start: b(k) = 3*a(k-1)*l(k-1) + b(k-1)
   tmp = zeros(1,4*nk); tmp(l) = [-3*lk(kk1) -1 0 0 0 1 0 0 0 0 0 0];
   A = [A ; tmp];
   B = [B ; 0];
  end
 end
end
C = reshape(A\B,4,nk).';

% ppform in linear coordinate
sk = [0;cumsum(lk)];
pp  = mkpp(sk,C);
ppd = mkpp(sk,[3*C(:,1) 2*C(:,2) C(:,3)]);

if nargin == 0
 plot(ppval(linspace(0,sk(end),101),pp),'.'), axis equal
end
