function [mat]=fantom_InNodes(vtx,simp,srf,mat_ref)
close all

% Set the background conductivity to be equal to the reference case
mat = mat_ref;
    % The range of base radius
        rad=[max(vtx(:,1))-min(vtx(:,1)),max(vtx(:,2))-min(vtx(:,2))];
    % The range of height
        height = [max(vtx(:,3))-min(vtx(:,3))];

        InhNmb = input('Please enter the number of inhomogeneities: ');

        for loop2 = 1:InhNmb

            %The Radius of the base of cylinder
            RInh(loop2) = input (['Radius of the base ([dx,dy]= [' num2str(rad(1)) ',' num2str(rad(2)) '] ): ']);

            %The Height of the cylinder
            HInh(loop2) = input (['Height of the cylinder ([dz]= [' num2str(height(1)) '] ):']);

            % Co-ordinates of the centre of the inhomogeneity
            CntInh(loop2,1) = input(['CentrumInh X co-ordinate (max=' num2str(max(vtx(:,1))) ', min=' num2str(min(vtx(:,1))) '): ']);
            CntInh(loop2,2) = input(['CentrumInh Y co-ordinate (max=' num2str(max(vtx(:,2))) ', min=' num2str(min(vtx(:,2))) '): ']);
            CntInh(loop2,3) = 0.0;


            % Co-ordinates of the top and bottom of the inhomogeneity
            TopInh(loop2) = input(['top Z co-ordinate (max=' num2str(max(vtx(:,3))) ', min=' num2str(min(vtx(:,3))) '): ']);
            BottomInh(loop2) = input(['bottom Z co-ordinate (max=' num2str(TopInh(loop2)) ', min=' 	num2str(min(vtx(:,3))) '): ']);

            % Conductivity of the inhomogenity
            CndInh(loop2) = input(['Conductivity within inhomogeneity : ']);

            % Initiate tetrahedra's counter
            SzInh(loop2) = 0;

            for loop1 = 1:size(vtx,1)

                % Return XYZ co-ordinates for the centroid of this tetrahedral element
                Cntrd(loop1,:) = vtx(loop1,:);

                % Calculate the distance between this centroid and the centre of the inhomogeneity
                DInhCntrd(loop2,loop1) = sqrt(sum((Cntrd(loop1,1:2)-CntInh(loop2,1:2)).^2));
              
                if ((Cntrd(loop1,3) <= TopInh(loop2)) && (Cntrd(loop1,3) >= BottomInh(loop2)))
               
                    if (DInhCntrd(loop2,loop1) <= RInh(loop2))

                        mat(loop1)=CndInh(loop2); % Change conductivity within inhomogeneity

                        SzInh(loop2) = SzInh(loop2)+1;

                    end

                end


            end

            disp(['Inhomogeneity (Cylinder) is comprised of ' num2str(SzInh(loop2)) ' tetrahedra'])
        end
