%GSZRC  Ancillary data for GSZR Poisson solver
% [CX,CQ,CR,CS,CI,CO] = GSZRC(R,DR,DZ,NR,NZ) returns ancilliary data to be
% used by GSZR,GSZRMEX. See also GSZRMEX.
%
% R: inner grid (ry)
% DR, DZ: r,z grid spacing.
% NR,NZ: inner (nry,nrz) grid dimensions.
%
% For details, see: [J-M.Moret et al. Fus.Eng.Des 2015] Section 4.2
%              and ﻿[F. Hofmann, Comput. Phys. Commun. 48(2) 1988].
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function [cx,cq,cr,cs,ci,co] = gszrc(r,dr,dz,nr,nz)
  assert(round(log2(nz+1)) == log2(nz+1),'log2(nz+1) must be an integer');
  assert(iscolumn(r),'r must be a column vector');
  
	c2 = (dz/dr)^2 * r;
	ce = c2 ./ (r + dr/2); co = ce(end)*0.5;
	cw = c2 ./ (r - dr/2); ci = cw(  1)*0.5;
	c2 = 2 + ce + cw;
	two = zeros(1,nz);
	t = sqrt(2);
	for k = 2.^(log2(nz+1)-2:-1:0)
		two(k) = t;
		l = 3*k:2*k:(nz+1)/2-1;
		two(l) = (two(l+k)+two(l-k))/t;
		t = sqrt(2+t);
	end
	two((nz+1)/2+1:end) = -fliplr(two(1:(nz+1)/2-1));
	cr = zeros(nr,nz);
	cq = c2-two;
	for i = nr-1:-1:1
		cr(i,:) = ce(i)./cq(i+1,:);
		cq(i,:) = cq(i,:) - cr(i,:).*cw(i+1);
	end
	cq = 1./cq;
	cs = cw.*cq;
	cx = pi*mu0*dz/dr*r;
end
