Moving to a unique version of LIUQE IDS to include all cases

IDS/gen_filament.m

-function  [discretizationr, discretizationz, Twc] = eval_trap(vvdata, num, type, varargin )
-% Evaluate the trapezoid from vv data structure, and computes the
-% discretization
-% type = 2 -> rectangle description
-% type = 3 -> oblique description
-% type = 4 -> contour of the coils are specified varargin = [r, z ]
-% contours
-% vv data contains R, Z, dR, dZ for type = 2.  
-% vv data contains R, Z, dR, dZ, AC, AC2 for type = 3 where AC, AC2 are the
-% angles of the parallelogram. Refer to GA document for specifications.
+function  [discretizationr, discretizationz, Twc] = gen_filament(vvdata, num, type_coil_description,varargin)
+% Generate filamentary(windings) discretization from different geometrical coil descriptions
 
-if mod(num,1)~=0, warning('The number of turns must be an integer. Ceil operation on it'); end
-num = ceil(num);
+% 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
+switch type_coil_description
     case 1 % outline option
-        disp('Not yet implemented. It will contain outline option IMAS')
-        return
+        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
+    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
-        disp('case not defined')
-        return
+        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 windings to active coil
-for ii = 1:numel(AC) % loop all the coils
+for ii = 1:numel(AC) % loop all the coils given vvstructure
+    % Find 4 extrema of the parallelogram
+    %   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 
         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 
         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']);
     
-    %fixed_area = 0.01; %Decide the area represented by each single coils
-    %area_of_element = abs(AB(1)*AD(2)-AB(2)*AD(1)); %Evaluate the closer area to the element
+    % 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);
     
-    %num = round(area_of_element/fixed_area); %Find the number of integer
-    %filaments closer to a fixed desired area for each element.
-    
-
-    if numel(num) == 1 %selfsimilar equally spaced
-        
-        pointx = zeros(num*num, 1);
-        pointy = zeros(num*num, 1);
-        
-        % Mode with selfsimilar shape
-        for kk=1:num
-            for jj=1:num
-                pointx((kk-1)*num+jj) = (2*kk-1)/(2*num); %uniform grid in the deformed space defined by the parallelogram
-                pointy((kk-1)*num+jj) = (2*jj-1)/(2*num);
-                
-                tmp = S*[pointx((kk-1)*(num)+jj), pointy((kk-1)*(num)+jj)]'; % evaluate component in the r,z system of reference
-                discretizationr(end+1) = pointr(ii,1)+ tmp(1) ; % translate to the original point
-                discretizationz(end+1) = pointz(ii,1) + tmp(2);
-            end
+    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
-        
-        % Creating the grouping matrix
-        tmp = zeros(num*num, numel(AC));
-        tmp(:,ii) = 1;
-        
-        Twc = [Twc ; tmp];
-        
-    else %Depending on the turns and same similar aspect ration
-        ratio = sqrt(dot(AB,AB)/dot(AD,AD));
-        
-        nx = ceil(sqrt(num(ii)/ratio));
-        ny = ceil(nx*ratio);
-        
-        if ny>=2
-            pointy = linspace(1/(ny-1)/2, 1- 1/(ny-1)/2, ny);
-        else
-            pointy = 0.5;
-        end
-        
-        if nx>=2
-            pointx = linspace(1/(nx-1)/2, 1- 1/(nx-1)/2, nx);
-        else
-            pointx = 0.5;
-        end
-        
-
-        [ pointxx, pointyy ] = meshgrid(pointx, pointy);
-        
-         for kk=1:ny
-            for jj=1:nx
-                
-                tmp = S*[pointxx(kk,jj); pointyy(kk,jj)]; % evaluate component in the r,z system of reference
-                discretizationr(end+1) = pointr(ii,1)+ tmp(1) ; % translate to the original point
-                discretizationz(end+1) = pointz(ii,1) + tmp(2);
-            end
-         end
-         
-         % Creating the grouping matrix
-         tmp = zeros(nx*ny, numel(num));
-         tmp(:,ii) = 1; 
-         Twc = [Twc ; tmp];
-         
     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 respect the input data.
+    Twc(end+1:end+num1*num2,end+1) = 1*sign(num(ii))/(num1*num2)*num(ii);
 end
 
-
-%% Plot the coils discretization from the original data
-%{
-figure
-hold on
-axis equal
-for ii = 1:numel(AC) % loop all the coils
-    patch(pointr(ii,:),pointz(ii,:),'red');
-end
-scatter(discretizationr, discretizationz, '.');
-%}
+% 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