function [r,z] = bouD(ntheta,R,Z,a,k1,k2,d1,d2,l1,l2,p1,p2)
% Analytical TCV (doublet) boundary parametrization with indentations, from 
% from Degtyarev EPS 1994 
% 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
% p1, p2: inner.outer indentation
%
% examples:
% % Simple shape with up-down asymmetric triangularity
%  [r,z] = bouD(31,0.85,0,0.25,1.4,1.4,0.5,0.1,0.1,0.1,0,0); plot(r,z,'o');
% % Doublet contour a la Degtyarev
%  [ra,za] = bouD(31,0.8,0,0.5,1,1,0,0,-0.1,-0.1,0.9,0.9); plot(ra,za,'o');
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

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

p = zeros(size(theta)); p(theta<pi) = p1; p(theta>pi)=p2;

r = R+a*cos(theta+d.*sin(theta)-l.*sin(2*theta)).*(1-p.*cos(theta).^2);
z = Z+a*k.*sin(theta);

end