Clean up versio of the function to generate the discretized filemants

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.
+
+if mod(num,1)~=0, warning('The number of turns must be an integer. Ceil operation on it'); end
+num = ceil(num);
+
+R = vvdata(1, :);
+Z = vvdata(2, :);
+dR = vvdata(3, :);
+dZ = vvdata(4, :);
+
+switch type
+    case 1 % outline option
+        disp('Not yet implemented. It will contain outline option IMAS')
+        return
+    case 2 % rectangle
+        AC = 0*ones(size(R));
+        AC2 = 90*ones(size(R));
+        
+    case 3 % oblique
+        AC = vvdata(5, :);
+        AC2 = vvdata(6, :);
+        
+    case 4 % contour of the polygon
+        
+    otherwise
+        disp('case not defined')
+        return
+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
+    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
+    
+    
+    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];
+    
+    
+    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
+    
+    %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
+        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
+    
+    
+    
+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, '.');
+%}
+discretizationr = reshape(discretizationr, numel(discretizationr),1);
+discretizationz = reshape(discretizationz, numel(discretizationz),1);
+
+
+end
+
+