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

% Set the background conductivity to be equal to the reference case
mat = mat_ref;

disp('The inhomogenity type');
disp('1.Sphere');
disp('2.Cylinder');
disp('3.Cylinder with coaxial hole');
disp('4.Fluid level');
disp('5.Hopper Flow');
disp('6.Fluid level with bubbles');
disp('7.Sphere inside a cylinder');
choice = input('Choose one : ');

switch(choice)
    % Sphere
    case 1
        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))) '): ']);

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

            % Initiate vertices's 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(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,:)-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 ' (Sphere) is comprised of ' num2str(SzInh(loop2)) ' tetrahedra'])
        end

    case 2

        % 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;

            % 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(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,:)-CntInh(loop2,:)).^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

        
                % Cylinder with coaxial hole
            case 3
        
                disp('OK!!!');
        
        
                % The range of volume 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))];
        
                %single inhomogenity
        
                InhNmb = input('Please enter the number of inhomogeneities: ');
        
                loop2 = 1;
        
                for loop2 = 1:InhNmb
        
                    %The Outer Radius of the base of cylinder with coaxial hole
        
                    ROutInh(loop2) = input (['Outer Radius of the base ([dx,dy]= [' num2str(rad(1)) ',' 	num2str(rad(2)) '] ): ']);
        
                    %The Inner Radius of the base of cylinder with coaxial hole
        
                    RInnInh(loop2) = input (['Inner Radius of the base ([dx,dy]= [' num2str(rad(1)) ',' 			num2str(rad(2)) '] ): ']);
        
                    %The Height of the cylinder with coaxial hole
        
                    HInh(loop2) = input (['Height of the cylinder ([dz]= [' num2str(height(1)) '] ):']);
        
        
                    % 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) = 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 inhomogeneity
        
                    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) <= ROutInh(loop2) && DInhCntrd(loop2,loop1) >= 		RInnInh(loop2))
        
                                mat(loop1)=CndInh(loop2); % Change conductivity within inhomogeneity
        
                                SzInh(loop2) = SzInh(loop2)+1;
        
                            end
        
                        end
        
        
                    end
        
                    disp(['Inhomogeneity (Cylinder with coaxial hole) is comprised of ' num2str(SzInh(loop2)) ' tetrahedra'])

                end
        
            case 4
                'Fluid Pipe'
                InhNmb = input('Please enter the number of inhomogeneities: ');
                for loop2 = 1:InhNmb
                    InhCnt = num2str(loop2);
                    FluidLevel = input (['Fluid Level ([0-100%]: ']);
                    Level=min(vtx(:,2))+(FluidLevel/100)*(max(vtx(:,2))-min(vtx(:,2)));
                    % 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))];
                    RInh(loop2) = rad(2)/2;
                    %The Height of the cylinder
                    HInh(loop2) = height(1);
                    % Co-ordinates of the centre of the inhomogeneity
                    CntInh(loop2,1) = 0.0;
                    CntInh(loop2,2) = 0.0;
                    CntInh(loop2,3) = 0.0;
                    % Co-ordinates of the top and bottom of the inhomogeneity
                    TopInh(loop2) = max(vtx(:,3));
                    BottomInh(loop2) = min(vtx(:,3));
                    % Conductivity of the inhomogenity
                    CndInh(loop2) = input(['Conductivity within inhomogeneity : ']);
                    % Initiate tetrahedra's 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(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
                        if Cntrd(loop1,2) < Level
                            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
        
        
                    end
        
                    disp(['Inhomogeneity is comprised of ' num2str(SzInh(loop2)) ' tetrahedra'])
                end
        
        
        %
        %     case 5
        %
        %         % The range of volume 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))];
        %         %single inhomogenity
        %         loop2 = 1;
        %         % 	%The Outer Radius of the base of cylinder with coaxial hole
        %         % 	ROutInh(loop2) = input (['Outer Radius of the base ([dx,dy]= [' num2str(rad(1)) ',' 	num2str(rad(2)) '] ): ']);
        %         % 	%The Inner Radius of the base of cylinder with coaxial hole
        %         % 	RInnInh(loop2) = input (['Inner Radius of the base ([dx,dy]= [' num2str(rad(1)) ',' 			num2str(rad(2)) '] ): ']);
        %         % 	%The Height of the cylinder with coaxial hole
        %         HInh(loop2) = input ([' Initial height of the hopper flow ([dz]= [' num2str(height(1)) '] ):']);
        %         % 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) = 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 inhomogeneity
        %         CndInh(loop2) = input(['Conductivity : ']);
        %         % Initiate tetrahedra's counter
        %         SzInhVol(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 mode
        %
        %         % HOPPER VOLUME
        %
        %
        %         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,:),:));
        %             min_x = min(vtx(:,1));
        %             max_x = max(vtx(:,1));
        %             min_y = min(vtx(:,2));
        %             max_y = max(vtx(:,2));
        %
        %             for x=min_x:max_x
        %                 for y=min_y:max_y
        %                     z_factor = - max_x/(x^2+y^2)
        %                     pause
        %                 end
        %             end
        %             % 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) <= ROutInh(loop2) && DInhCntrd(loop2,loop1) >= RInnInh(loop2))
        %             %                 mat(loop1)=CndInh(loop2); % Change conductivity within inhomogeneity
        %             %                 SzInh(loop2) = SzInh(loop2)+1;
        %             %             end
        %             %         end
        %         end
        %
        %         disp(['Inhomogeneity (Hopper Flow Model) is comprised of ' num2str(SzInh(loop2)) ' tetrahedra'])
        %         figure
        %         trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3))
        %         axis image
        %         view(54,70)
        %         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 gouraud
        %         pause(3)
        %
            case 6
                'Fluid Pipe with Bubbles'
                %
                for loop2 = 1:1
                    InhCnt = num2str(loop2);
                    FluidLevel = input (['Fluid Level ([0-100%]: ']);
                    Level=min(vtx(:,2))+(FluidLevel/100)*(max(vtx(:,2))-min(vtx(:,2)));
                    % 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))];
                    RInh(loop2) = rad(2)/2;
                    %The Height of the cylinder
                    HInh(loop2) = height(1);
                    % Co-ordinates of the centre of the inhomogeneity
                    CntInh(loop2,1) = 0.0;
                    CntInh(loop2,2) = 0.0;
                    CntInh(loop2,3) = 0.0;
                    % Co-ordinates of the top and bottom of the inhomogeneity
                    TopInh(loop2) = max(vtx(:,3));
                    BottomInh(loop2) = min(vtx(:,3));
                    % Conductivity of the inhomogenity
                    CndInh(loop2) = input(['Conductivity within inhomogeneity : ']);
                    % Initiate tetrahedra's 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(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
                        if Cntrd(loop1,2) < Level
                            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
        
        
                    end
        
                    rad=[max(vtx(:,1))-min(vtx(:,1)),max(vtx(:,2))-min(vtx(:,2)),max(vtx(:,3))-min(vtx(:,3))];
                    InhSphNmb = input('Please enter the number of inhomogeneities: ');
                    for loop2 = 1:InhSphNmb
                        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))) '): ']);
                        % Conductivity of the inhomogeneity
                        CndInh(loop2) = input(['Conductivity within inhomogeneity ' InhCnt ': ']);
                        % Initiate tetrahedra's 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(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,:)-CntInh(loop2,:)).^2));
                            if DInhCntrd(loop2,loop1)<=RInh(loop2)
                                mat(loop1)=CndInh(loop2); % Change conductivity within inhomogeneity
                                SzInh(loop2) = SzInh(loop2)+1;
                            end
                        end
                    end  
        
                    disp(['Inhomogeneity is comprised of ' num2str(SzInh(loop2)) ' tetrahedra'])
                end
                
             case 7   
                
                 disp('OK!!!');
        
        
                % The range of volume 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))];
        
                %single inhomogenity
        
                InhNmb = input('Please enter the number of inhomogeneities: ');
        
                loop2 = 1;
        
                for loop2 = 1:InhNmb
        
                    %The Outer Radius of the base of cylinder with coaxial hole
        
                    ROutInh(loop2) = input (['Outer Radius of the base ([dx,dy]= [' num2str(rad(1)) ',' 	num2str(rad(2)) '] ): ']);
        
                    %The Inner Radius of the base of cylinder with coaxial hole
        
                    RInnInh(loop2) = input (['Inner Radius of the base ([dx,dy]= [' num2str(rad(1)) ',' 			num2str(rad(2)) '] ): ']);
        
                    %The Height of the cylinder with coaxial hole
        
                    HInh(loop2) = input (['Height of the cylinder ([dz]= [' num2str(height(1)) '] ):']);
        
        
                    % 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) = 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 inhomogeneity
        
                    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) <= ROutInh(loop2) && DInhCntrd(loop2,loop1) >= 		RInnInh(loop2))
        
                                mat(loop1)=CndInh(loop2); % Change conductivity within inhomogeneity
        
                                SzInh(loop2) = SzInh(loop2)+1;
        
                            end
        
                        end
        
        
                    end
        
                    disp(['Inhomogeneity (Cylinder with coaxial hole) is comprised of ' num2str(SzInh(loop2)) ' tetrahedra'])

                end
                
                    rad=[max(vtx(:,1))-min(vtx(:,1)),max(vtx(:,2))-min(vtx(:,2)),max(vtx(:,3))-min(vtx(:,3))];
                    InhSphNmb = input('Please enter the number of inhomogeneities: ');
                    for loop2 = 1:InhSphNmb
                        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))) '): ']);
                        % Conductivity of the inhomogeneity
                        CndInh(loop2) = input(['Conductivity within inhomogeneity ' InhCnt ': ']);
                        % Initiate tetrahedra's 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(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,:)-CntInh(loop2,:)).^2));
                            if DInhCntrd(loop2,loop1)<=RInh(loop2)
                                mat(loop1)=CndInh(loop2); % Change conductivity within inhomogeneity
                                SzInh(loop2) = SzInh(loop2)+1;
                            end
                        end
                    end  
        
                    disp(['Inhomogeneity is comprised of ' num2str(SzInh(loop2)) ' tetrahedra'])

    otherwise
        disp('Bad choice :(');

end