function [mat]=ECT3D_FEMCreate_Inhomogenity_Cylinder(vtx,simp,srf,mat_ref,from_min,to_max);
%[mat]=ECT3D_FEMCreate_Inhomogenity_Cylinder(vtx,simp,srf,mat_ref,from_min,to_max);
%
%This function create ingomegenity inside object.
%mat - new material with the inhomogenity
%
%mat_ref - base material
%from_min,to_max - only those elements inside this range can be changed

close all


if nargin <6
    from_min=1;
    to_max=max(size(mat_ref));
end
% Set the background conductivity to be equal to the reference case
mat = mat_ref;

rad=[max(vtx(:,1))-min(vtx(:,1)),max(vtx(:,2))-min(vtx(:,2)),max(vtx(:,3))-min(vtx(:,3))];

InhNmb = input('Please enter the number of inhomogeneities: ');
for loop2 = 1:InhNmb
    InhCnt = num2str(loop2);
    RInh(loop2) = input (['Radius of inhomogeneity ' num2str(loop2) ' ([dx,dy,dz]= [' num2str(rad(1)) ',' num2str(rad(2)) ',' num2str(rad(3)) '] ): ']);
    % Co-ordinates of the centre of the inhomogeneity
    CntInh(loop2,1) = input(['X co-ordinate (max=' num2str(max(vtx(:,1))) ', min=' num2str(min(vtx(:,1))) '): ']);
    CntInh(loop2,2) = input(['Y co-ordinate (max=' num2str(max(vtx(:,2))) ', min=' num2str(min(vtx(:,2))) '): ']);
    %CntInh(loop2,3) = input(['Z co-ordinate (max=' num2str(max(vtx(:,3))) ', min=' num2str(min(vtx(:,3))) '): ']);
    CndInh(loop2) = input(['Conductivity within inhomogeneity ' InhCnt ': ']);
    % Initiate tetradra counter
    SzInh(loop2) = 0;
    % Create a matrix Cntrd comprised of the XYZ co-ordinates, with
    % respect to the origin, of the centroid of each tetrahedral
    % element of the wire mesh model
    for loop1 = 1:size(simp,1)
        if loop1<from_min | loop1>to_max 
            continue;
        end
        % Return XYZ co-ordinates for the centroid of this tetrahedral element
        Cntrd(loop1,:) = mean(vtx(simp(loop1,:),1:2));
        % Calculate the distance between this centroid and the centre of the inhomogeneity
        DInhCntrd(loop2,loop1) = sqrt(sum((Cntrd(loop1,:)-CntInh(loop2,:)).^2));
        if DInhCntrd(loop2,loop1)<=RInh(loop2)
            mat(loop1)=CndInh(loop2); % Change conductivity within inhomogeneity
            SzInh(loop2) = SzInh(loop2)+1;
        end
    end
    disp(['Inhomogeneity ' InhCnt ' is comprised of ' num2str(SzInh(loop2)) ' tetrahedra'])
end
DInhCntrd = DInhCntrd';

% % Plot the inhomogeneities within the uniform bulk material
% figure
% trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3))
% axis image
% view(45,10)
% set(gcf,'Name','Inhomogeneities in the Bulk','Colormap',[0,0,0])
% hidden off
% hold on
% ECT3D_FEMRepaint_Inho(mat,mat_ref,vtx,simp)
% camlight lefth
% lighting flat
% pause(3)