function [MasterC1,masterC_list1,weightC1,ia1] = SelectContact(inp,g1,r,numpos)

%% SELECTION OF THE NODES FOR CONTACT STIFFNESS AT EACH ANGULAR POSITION AND TOOTH
%  Plot di controllo selezione nodi
    figure, hold on;
    plot3(g1.FEM.nodes(:,2),g1.FEM.nodes(:,3),g1.FEM.nodes(:,4),'.k');
    xlabel('X [mm]');
ylabel('Y [mm]');
zlabel('Z [mm]');
    axis equal,
p1=cell(numpos,inp.params.sel_teeth); %initialize as empty cell array 
for iii= 1:numpos % loop for the angular position counter
for tk=2:4 %loop for all teeth - inp.params.sel_teeth
    for lk=2:inp.params.slices-1 %loop for alls slices
        % First and last slice are excluded otherwise i get points on the
        % face of the gear - To be addressed in the future
        % Pinion nodes (index is lk-1 because i'm starting from 2)
if not(isnan(r(1).gearpair.CNall{1,iii}{iii,tk}(lk,2))) && ...
        r(1).gearpair.CNall{1,iii}{iii,tk}(lk,2)<=size(g1.rigid_contact.xmat,1) &&...
        r(1).gearpair.CNall{1,iii}{iii,tk}(lk,4)<=inp.params.contact_tol_archi*180/pi

        p1{iii,tk}(lk-1,1) =  g1.rigid_contact.xmat(r(1).gearpair.CNall{1,iii}{iii,tk}(lk,g1.number+1),r(1).gearpair.CNall{1,iii}{iii,tk}(lk,1),tk);
        p1{iii,tk}(lk-1,2) =  g1.rigid_contact.ymat(r(1).gearpair.CNall{1,iii}{iii,tk}(lk,g1.number+1),r(1).gearpair.CNall{1,iii}{iii,tk}(lk,1),tk);
        p1{iii,tk}(lk-1,3) =  g1.rigid_contact.zmat(r(1).gearpair.CNall{1,iii}{iii,tk}(lk,g1.number+1),r(1).gearpair.CNall{1,iii}{iii,tk}(lk,1),tk);
end     
end %loop for all slices
end %loop for all teeth

for tk=1:inp.params.sel_teeth %loop for all teeth - inp.params.sel_teeth
    if not(isempty(p1{iii,tk}))
    [nod1{iii,tk},distC1{iii,tk}]=contact2fem([g1.FEM.surfnodes(:,2),g1.FEM.surfnodes(:,3),g1.FEM.surfnodes(:,4)],...
        p1{iii,tk});
    [masterC_list1{iii,tk},ia1{iii,tk},~]=unique(nod1{iii,tk},'stable');

    % ensure that masterC_list is a column (when there are few nodes unique
    % sometimes does not perform reshaping)
    if size(masterC_list1{iii,tk},2)>1
        tempa=masterC_list1{iii,tk};
        masterC_list1{iii,tk}=reshape(masterC_list1{iii,tk},[],1);
    end

    % Calculation of the weight to be assigned to each node to spread the
    % contact stiffness over the contact area
    for sl=1:size(distC1{iii,tk},1)
        di=distC1{iii,tk}(sl,:);
        w=di/sum(di);
        weightC1{iii,tk}(sl,:)=w;
    end
    
   
    plot3(g1.FEM.surfnodes(masterC_list1{iii,tk},2),g1.FEM.surfnodes(masterC_list1{iii,tk},3),g1.FEM.surfnodes(masterC_list1{iii,tk},4),'xr'); hold on,
    plot3(p1{iii,tk}(:,1),p1{iii,tk}(:,2),p1{iii,tk}(:,3),'ob');
    end
end % end loop for all teeth
end
% All unique nodes needed to attach the contact stiffness, in ascending
% order. masterC_list will be used for each angular position since it's
% divided by angular position and tooth number
tempc=reshape(masterC_list1,[],1);
tempc=tempc(~cellfun('isempty',tempc));
MasterC1=unique(cell2mat(tempc));
