function [r,z] = bouL(ntheta,R,Z,a,k1,k2,d1,d2,l1,l2)
% Lemniscate boundary
% R,  Z : center R,Z
% a     : minor radius (overall size)
% k1, k2: top/bottom elongation
% d1, d2: top/bottom triangularity
% l1, l2: top/bottom squareness
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

th = linspace(0,2*pi,ntheta+1)'; th = th(1:ntheta);
k = 0.5*(k1 + k2 + (k1-k2)*sin(th));
d = 0.5*(d1 + d2 + (d1-d2)*sin(th));
l = 0.5*(l1 + l2 + (l1-l2)*sin(th));

thh = th-pi/2;
r = R+(2*a.*sin(thh).*cos(thh))./(1+sin(thh).^2);
z = Z+(2*a.*k.*cos(thh+d.*sin(th)+l.*sin(2*th)))./(1+sin(thh).^2);

end