function [zA,rA,FA,dz2FA,dr2FA,drzFA,zX,rX,FX,dz2FX,dr2FX,drzFX] = ...
  asxycs(Fx,zA,rA,zX,rX,L)
% Refine asxy result using cubic spline interpolation
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

npar = L.csnthreads;
rk = L.taur;
zk = L.tauz;

% spline coefficients
c = L.Mr*(L.Mz*Fx).';

nA = numel(rA);
rs = [rA;rX]; zs = [zA;zX]; ns = numel(rs);
Fs = zeros(ns,1); dr2Fs = Fs; dz2Fs = Fs; drzFs = Fs;

niter = L.P.itercs;
tol = L.P.tolcs;

% Gauss-Newton iteration
for kiter = 1:niter
  [~,gr,gz,hrr,hrz,hzz] = bsp2deval(rk,zk,c,rs,zs,[1,2],npar);
  idh = 1./(hrr.*hzz-hrz.^2);
  er = (hzz.*gr-hrz.*gz).*idh;
  ez = (hrr.*gz-hrz.*gr).*idh;
  rs = rs - er; zs = zs - ez;
  if max(max(abs(er)),max(abs(ez))) < tol
    [Fs,~,~,dr2Fs,drzFs,dz2Fs] = bsp2deval(rk,zk,c,rs,zs,[0,2],npar);
    break
  end
end

% Dispatch
rA    =    rs(1:nA,1);
zA    =    zs(1:nA,1);
FA    =    Fs(1:nA,1);
dr2FA = dr2Fs(1:nA,1);
drzFA = drzFs(1:nA,1);
dz2FA = dz2Fs(1:nA,1);
rX    =    rs(nA+1:ns,1);
zX    =    zs(nA+1:ns,1);
FX    =    Fs(nA+1:ns,1);
dr2FX = dr2Fs(nA+1:ns,1);
drzFX = drzFs(nA+1:ns,1);
dz2FX = dz2Fs(nA+1:ns,1);

end
