function  [discretizationr, discretizationz, Twc] = gen_filament(vvdata, num, type_coil_description,varargin)
% Generate filamentary(windings) discretization, and grouping matrix, from different geometrical coil descriptions

% type_coil_description = 2 -> rectangle description
% type_coil_description = 3 -> oblique description
% type_coil_description = 4 -> contour of the coils are specified varargin = [r, z ] contours
% vvdata =  R, Z, dR, dZ for type_coil_description = 2.  
% vvdata =  R, Z, dR, dZ, AC, AC2 for type_coil_description = 3 where AC, AC2 are the angles of the parallelogram. 

R = vvdata(1, :);
Z = vvdata(2, :);
dR = vvdata(3, :);
dZ = vvdata(4, :);
switch type_coil_description
    case 1 % outline option
        error('Filamentary discretization for ''outline'' coils description in ids type_coil_description not yet implmented')
    case 2 % rectangle
        AC = 0*ones(size(R));
        AC2 = 90*ones(size(R));
    case 3 % oblique (parallelogram)
        AC = vvdata(5, :);
        AC2 = vvdata(6, :);
    case 4 % contour of the polygon
        error('Filamentary discretizaion for ''polygon'' coils description in  ids typ not yet implemented')
    otherwise
        error('Coils description not supported for filamenteray discretization ')
end
        
% Initialize the extrema of the trapezoid
pointr = zeros(numel(AC),4);
pointz = zeros(numel(AC),4);
discretizationr = [];
discretizationz = [];
Twc = []; % Grouping matrix from coils to windings Iw = Twc*Ic
for ii = 1:numel(AC) % loop all the coils given vvstructure
    % Find 4 extrema of the parallelogram. Working both for rectangular description and parallelogram description
    %   D----C
    %  /    /
    % A----B
    if  AC2(ii) ~= 0 && AC(ii) ==0 
        DH = dZ(ii)/tand(AC2(ii));
        pointr(ii,1) = R(ii) -DH/2 -dR(ii)/2;
        pointr(ii,2) = R(ii) +DH/2 -dR(ii)/2;
        pointr(ii,3) = R(ii) +DH/2 +dR(ii)/2;
        pointr(ii,4) = R(ii) -DH/2 +dR(ii)/2;
        pointz(ii,1) = Z(ii) -dZ(ii)/2;
        pointz(ii,2) = Z(ii) +dZ(ii)/2;
        pointz(ii,3) = Z(ii) +dZ(ii)/2;
        pointz(ii,4) = Z(ii) -dZ(ii)/2;
    elseif AC(ii) ~= 0 && AC2(ii) ==0 % Parallelogramm
        DH = dR(ii)*tand(AC(ii));
        pointr(ii,1) = R(ii) - dR(ii)/2;
        pointr(ii,2) = R(ii) - dR(ii)/2;
        pointr(ii,3) = R(ii) + dR(ii)/2;
        pointr(ii,4) = R(ii) + dR(ii)/2;
        pointz(ii,1) = Z(ii) - DH/2- dZ(ii)/2;
        pointz(ii,2) = Z(ii) - DH/2 +dZ(ii)/2;
        pointz(ii,3) = Z(ii) + DH/2 +dZ(ii)/2;
        pointz(ii,4) = Z(ii) + DH/2 -dZ(ii)/2;
    elseif AC2(ii) == 0 && AC(ii) ==0 % Rectangular description 
        % 2----3
        % |    |
        % 1----4
        pointr(ii,1) = R(ii)  -dR(ii)/2;
        pointr(ii,2) = R(ii)  -dR(ii)/2;
        pointr(ii,3) = R(ii)  +dR(ii)/2;
        pointr(ii,4) = R(ii)  +dR(ii)/2;
        pointz(ii,1) = Z(ii) -dZ(ii)/2;
        pointz(ii,2) = Z(ii) +dZ(ii)/2;
        pointz(ii,3) = Z(ii) +dZ(ii)/2;
        pointz(ii,4) = Z(ii) -dZ(ii)/2;
    end
    % Define base vector of the 2D deformed space for the parallelogram
    AD = [pointr(ii,4)-pointr(ii,1), pointz(ii,4)-pointz(ii,1)];
    AB = [pointr(ii,2)-pointr(ii,1), pointz(ii,2)-pointz(ii,1)];
    ex = [1,0];
    ey = [0,1];
    
    % Rotation matrix containing the scalar product of the deformed space
    % and orthogonal space
    S = ([ex*AD', ex*AB'; ey*AD', ey*AB']);
    
    % Try to use a number of filaments that respect the shape of the parallelogram
    ratio = sqrt((AD*AD'))/sqrt((AB*AB'));
    num1 = ceil(sqrt(num(ii)*ratio));
    num2 = ceil(num1/ratio);
    
    pointx = zeros(num1*num2, 1);
    pointy = zeros(num1*num2, 1);
    for kk=1:num1
        for jj=1:num2
            pointx((kk-1)*num1+jj) = (2*kk-1)/(2*num1); %uniform grid in the deformed space defined by the parallelogram
            pointy((kk-1)*num2+jj) = (2*jj-1)/(2*num2);
            tmp = S*[pointx((kk-1)*(num1)+jj), pointy((kk-1)*(num2)+jj)]'; % evaluate component in the r,z system of reference
            discretizationr(end+1) = pointr(ii,1) + tmp(1) ; % translate to the original central location
            discretizationz(end+1) = pointz(ii,1) + tmp(2);  % translate to the original central location
        end
    end
    % Only the number of turns is specified but not the position of
    % each individual winding. The function eval_trap generates a winding
    % distrubution that nicely resamble the geometry of the coil section,
    % but will not necessarily generate a number of windings = number of
    % turns. The following normalization garanties that the total current
    % flowing in the cross section of coil Iw*nw = Ia*nturns respects the input data.
    Twc(end+1:end+num1*num2,end+1) = 1*sign(num(ii))/(num1*num2)*num(ii);
end

% Reshape filaments(windings) into unique vector
discretizationr = reshape(discretizationr, numel(discretizationr),1);
discretizationz = reshape(discretizationz, numel(discretizationz),1);

% Plot the coils discretization
if ~isempty(varargin) && varargin{1} == 1
    figure
    hold on
    axis equal
    for ii = 1:numel(AC) % loop all the coils
        patch(pointr(ii,:),pointz(ii,:),'red');
    end
    scatter(discretizationr, discretizationz, '.');
end

end