classdef (Abstract) Shape3D<shape.Shape
    %SHAPE3D is the super class for  3d shapes
    %

    % Copyright 2022-2023 The MathWorks, Inc. 
    properties
       Dielectric = 'Air'
    end
    properties
         Tetrahedra
    end
    properties(Dependent)
        % Tetrahedra that forms the volume of shape
        TetrahedraShape
        MetalShape
    end
    properties(Access=protected)
       TetStruct = struct('TetrahedraShape',[],...
           'MetalShape',[]);
    end

    methods
        function obj = Shape3D()
            obj = obj@shape.Shape();
            obj.ShapeDimension='3D';
        end
    end

    methods
        function set.TetrahedraShape(obj,propVal)
              obj.TetStruct.TetrahedraShape = propVal;
        end
        function set.MetalShape(obj,propVal)
              obj.TetStruct.MetalShape = propVal;
        end
        function propVal = get.TetrahedraShape(obj)
              propVal = obj.TetStruct.TetrahedraShape ;
        end
        function propVal = get.MetalShape(obj)
              propVal = obj.TetStruct.MetalShape ;
        end
         function  set.Dielectric(obj,propVal)
            if isa(obj,'shape.Shape2D')
                if isequal(obj.Dielectric,propVal)
                    return;
                end
            end
            % validatestring(propVal,{'Air','FR4','Teflon','Foam','Polystyrene','Plexiglas',...
            %     'Fused quartz','E glass','RO4725JXR','RO4730JXR','TMM3',...
            %     'TMM4','TMM6','TMM10','TMM10i','Taconic RF-35'},class(obj),'Dielectric');
            met = DielectricCatalog;
            if ~strcmpi(propVal,'multidiel')
                validatestring(propVal,cellfun(@(x) (x), met.Materials.Name',...
                    'UniformOutput',false),class(obj),'Dielectric');
            end
            if isPropertyChanged(obj,obj.Dielectric,propVal)
                obj.Dielectric = propVal;
            end
        end
    end

    methods(Hidden)
        
        function [edgeLength,growthRate] = calculateMeshParams(obj,varargin)
            % Calculates the edge length required for meshing the shape


            if nargin>1
                if isprop(obj,'FeedWidth')
                    edgeLength1 = obj.FeedWidth*5;
                else
                    edgeLength1 = [];
                end
                edgeLength2 = varargin{1}/5;
                edgeLength = max([edgeLength1,edgeLength2]);
                growthRate = 1.75;
            else
                edgeLength       = [];
                growthRate = 1.75;
            end
        end

        function meshGenerator(obj,varargin)
            % Check if structure has changed and regenerate mesh if needed
            % Input is the shape(obj)
            % varargin include
            % edgeLength,growthRate,minEdgeLength,targetEdgeLength,
            % minFeatureEdgeLength and will be used as inputs for remesher
            
            if getHasStructureChanged(obj)||~isempty(obj.MesherStruct.Mesh.p)
                createGeometry(obj);
            end
            
            % Stores domain number before remeshing
            if ~isempty(obj.MesherStruct.Mesh.p)
                obj.MultiShapeID = obj.BooleanMultiShapeID;
            end

            % Mesh structure
            if numel(varargin)>1
                edgeLength  = varargin{1};
                growthRate  = varargin{2};
                setMeshEdgeLength(obj, edgeLength);
                setMeshGrowthRate(obj, growthRate);
                if numel(varargin) >2
                    minEdgeLength = varargin{3};
                    targetEdgeLength = varargin{4};
                    minFeatureEdgeLength = varargin{5};
                    setMeshMinContourEdgeLength(obj,minEdgeLength);
                    setTargetInteriorEdgeLength(obj,targetEdgeLength);
                    setSubDomainEdgeLength(obj,minFeatureEdgeLength);
                end
            else
                edgeLength  = getMeshEdgeLength(obj);
                growthRate  = getMeshGrowthRate(obj);
                minEdgeLength = getMinContourEdgeLength(obj);
                targetEdgeLength = 0;
                minFeatureEdgeLength=0;
            end
            % Assign points and triangles
            pr=obj.Points;tr=obj.Triangles;
            warnState = warning('off','MATLAB:triangulation:PtsNotInTriWarnId');
            if size(tr,2)==4
                TR = triangulation(tr(:,1:4),pr);
            elseif ~isempty(obj.MultiShapeID)
                TR = triangulation([tr(:,1:3) obj.MultiShapeID],pr);
            else
                TR = triangulation(tr(:,1:3),pr);
            end
            warning(warnState);
            p = TR.Points;
            t = TR.ConnectivityList;
            if size(t,2)==3
                %                 t(:,4) = 0;
                % For a normal box without operation,delete faces requires
                % layer ids to be present after remeshing,so that
                % corresponding faces can be removed.
                if isa(obj,'shape.Box')
                    t(:,4) =obj.InitialFaceID;
                else
                    t(:,4)=0;
                end
            end

            % Check mesh for triangles with small areas
            A = em.internal.meshprinting.meshareas(p, t);
            % Find triangles with areas approaching machine epsilon and remove
            Epsilon = 1e-12; %sqrt(eps);
            id = A<Epsilon;
            t(id,:) = [];
            FeedsAndVias=[];
            [ps,es,ts] = em.MeshGeometry.remesher(p,t,edgeLength,growthRate,minEdgeLength,...
                targetEdgeLength,minFeatureEdgeLength,  ...
                FeedsAndVias,[]);
            if isa(obj,'shape.Box')
                ts(4,end)=8;
                [ps,ts]=remapdomain(obj,p,t,ps',ts',FeedsAndVias);
            else
                [ps,ts]=remapdomain(obj,p,t,ps',ts',FeedsAndVias);
            end
            % Assign p,t from remesher
            obj.MesherStruct.Mesh.p=ps;
            obj.MesherStruct.Mesh.t=ts;
            obj.Edges = es;
            if isempty(obj.Tetrahedra)
                T           = [];
                epsr    = [];
                LTang = [];
            else
                % Only metal remeshes when shape has substrate and metal
                error(message('antenna:shapeerrors:NoRemeshForShapeWithSub'));
            end
            Mesh        = em.internal.makeMeshStructure(ps,ts,T,epsr,...
                LTang);
            % Save Mesh
            saveMesh(obj,Mesh);
        end

    end
    methods

         function   totobj = addSubstrate(tshape,tetshape)

             % This function adds dielectric behaviour to shape passed as
             % second argument. The first shape is treated as metal
             %
             % Example 1
             %
             % % Assign Substrate to a shape
             %
             % substratebox = shape.Box(Length = 0.03,Width = 0.025,Height = 1e-3,Center = [0 0 0.5e-3],Dielectric = 'Teflon',Color = [0.1000 0.7000 0.5000]);
             % Top = shape.Rectangle(Length = 0.0115,Width = 0.0130,Color = [0.8745 0.7255 0.2275]);
             % [~] = translate(Top,[0,0,1e-3]);
             % Bottom = shape.Rectangle(Length = 0.03,Width = 0.025,Color = [0.8745 0.7255 0.2275]);
             % AddTop = Bottom + Top;
             % ant = addSubstrate(AddTop,substratebox);
             % show(ant);
             
            % Parse Inputs
%             parserObj=inputParser;
%             addParameter(parserObj,'EpsilonR',1);
%             addParameter(parserObj,'LossTangent',0);
%             parse(parserObj,varargin{:});
%             epsr  = parserObj.Results.EpsilonR;
%             LossTg  = parserObj.Results.LossTangent;
%             tshape.EpsilonR=epsr;
%             tshape.LossTangent=LossTg;
            % tshape.Material.EpsilonR = tetshape.Material.EpsilonR;
            % tshape.Material.LossTangent = tetshape.Material.LossTangent;
            tetshape=copy(tetshape);
            % make tetrahedra
            % idx1 = (unique(tetshape.MultiShapeID));
            % idx2 = numel(unique(tshape.MultiShapeID));
            % idx1=idx1+idx2;
            totobj = xor(tshape,tetshape);
            % tt = [totobj.Triangles totobj.MultiShapeID];in_tri_idx = zeros(size(tt,1),1);
            % for n=1:numel(idx1)
            %     in_tri_idx(tt(:,4)==idx1(n))=1;
            % end
            try
                % Make tetra mesh
                [po, T,~,~,~,~] = em.MeshGeometry.tetramesher(totobj.Points',totobj.Triangles',MeshSize=0.1,GrowthRate=1.75);
            catch
                error(message("antenna:shapeerrors:TetMeshErrorFromShape"));

            end
            dTet=triangulation(T',po');
            % tshape.Tetrahedra=dTet;
            totobj.Tetrahedra=dTet;
            totobj.TetrahedraShape = copy(tetshape);
            totobj.MetalShape = copy(tshape);
            if isempty(tshape.MultiShapeID)
                totobj.MetalShape.MultiShapeID = ones(size(totobj.MetalShape.Triangles,1),1);
            end
            if isempty(tetshape.MultiShapeID)
                % totobj.TetrahedraShape.MultiShapeID = ones(size(totobj.TetrahedraShape.Triangles,1),1);
                 totobj.TetrahedraShape.MultiShapeID = ones(size(totobj.TetrahedraShape.Triangles,1),1)*...
                    (max(unique(totobj.MetalShape.MultiShapeID))+1);
            end
        end

        function varargout = rotate(obj,ang,axispt1,axispt2,varargin)
            % ROTATE function rotates the given shape by a certain
            % angle about the specified axis vector defined with
            % the two end points.
            %
            % % Examples:
            %
            % % Example 1:
            % % Create a shape and rotate it about the axis [0 0 0] and [0 1
            % % 0] by an angle of 90 degrees.
            % sh = shape.Box;
            % rotate(sh,90,[0 0 0],[0 1 0])
            %
            % % Example 2:
            % % Create a shape and rotate it by 45 degrees about the axis
            % % defined by the end points [ 0 0 0] and [0 0 1] and plot it on
            % % a figure handle
            %
            % sh = shape.Box;
            % figHandle = figure;
            % rotate(sh,45,[0 0 0],[0 0 1],'Parent',figHandle);

            % Inputs are shape(obj),angle(ang),axes
            % defined(axispt1,axispt2)
            % varargin include name value pair(Parent) which accepts a axis
            % or figure handle and output is plotted on that corresponding
            % axes or figure.
            
            

            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end
            if nargin==6
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{:});
                ParentObj=parserObj.Results.Parent;
                checkParentInput(obj,ParentObj);
            else
                ParentObj=[];
            end

            if getHasStructureChanged(obj)
                createGeometry(obj);
            end            
            % Rotate the object in 3D plane
            rotate3D(obj,ang,axispt1,axispt2);
            % Show final structure
            if nargout==0
                if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,'Parent',ParentObj);
                end
            else
                varargout{1} = obj;
            end
        end

        function varargout = rotateX(obj,ang,varargin)
            % rotateX rotates the shape along x-axis by specified angle
            %
            % % Examples:
            %
            % % Example1:
            % % Create a box shape and rotate it by 45
            % % degrees about the x- axis.
            %
            % r  = shape.Box;
            % rotateX(r,45);
            
            % Input is the 3d shape(obj)
            % varargin include name value pair (Parent) that accepts figure
            % or axes handle 


            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end

            if nargin==4
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{:});
                ParentObj=parserObj.Results.Parent;
                checkParentInput(obj,ParentObj);
            else
                ParentObj=[];
            end

            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            rotate3D(obj,ang,[0, 0, 0],[1, 0, 0]);

            % Show final structure
            if nargout==0
                if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,"Parent",ParentObj);
                end
            else
                varargout{1} = obj;
            end

        end

        function varargout = rotateY(obj,ang,varargin)
            % rotateY rotates the shape along y-axis by specified angle
            %
            % % Examples:
            %
            % % Example1:
            % % Create a box shape and rotate it by 45
            % % degrees about the y- axis.
            %
            % r  = shape.Box;
            % rotateY(r,45);
            
            % Input is the 3d shape(obj)
            % varargin include name value pair (Parent) that accepts figure
            % or axes handle 


            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end

            if nargin==4
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{:});
                ParentObj=parserObj.Results.Parent;
                checkParentInput(obj,ParentObj);
            else
                ParentObj=[];
            end

            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            rotate3D(obj,ang,[0 0 0],[0 1 0]);
            % Show final structure
            if nargout==0
                if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,"Parent",ParentObj);
                end
            else
                varargout{1} = obj;
            end
            
        end

        function varargout = rotateZ(obj,ang,varargin)
            % rotateX rotates the shape along z-axis by specified angle
            %
            % % Examples:
            %
            % % Example1:
            % % Create a box shape and rotate it by 45
            % % degrees about the z- axis.
            %
            % r  = shape.Box;
            % rotateZ(r,45);
            
            % Input is the 3d shape(obj)
            % varargin include name value pair (Parent) that accepts figure
            % or axes handle 


            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end
            if nargin==4
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{:});
                ParentObj=parserObj.Results.Parent;
                checkParentInput(obj,ParentObj);
            else
                ParentObj=[];
            end


            obj = rotate(obj,ang,[0, 0, 0],[0, 0, 1]);
            % Show final structure
            if nargout==0
                if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,"Parent",ParentObj);
                end
            else
                varargout{1} = obj;
            end
        end
        function varargout = translate(obj,vec,varargin)
            % translate function moves the shape to a new location 
            %
            % % Examples:
            %
            % % Example1:
            % % Create a box shape and translate it by it by [2 2 2]
            %
            % r  = shape.Box;
            % translate(r,[2 2 2]);
            
            % Input is the 3d shape(obj),translation vector(vec)
            % varargin include name value pair (Parent) that accepts figure
            % or axes handle 

            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end
            validateattributes(vec,{'numeric'},...
                    {'nonempty','nonnan','finite','real',...
                    'numel',3},'translate','Translation Vector',2);
            if ~isrow(vec)
                vec = vec';
            end
            if nargin==4
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{:});
                ParentObj=parserObj.Results.Parent;
                checkParentInput(obj,ParentObj);
            else
                ParentObj=[];
            end


            if getHasStructureChanged(obj)
                %             if isempty(obj.Points)
                createGeometry(obj);
            end
            translateOutOfXYPlane(obj,vec);
            if isprop(obj,'Center')
                %                 vec=orientGeom(obj,vec');
                obj.Center=obj.Center+vec;
                resetHasStructureChanged(obj);
            end
            % Show final structure
            if nargout==0
                if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,'Parent',ParentObj)
                end
            else
                varargout{1} = obj;
            end
        end

        function show(obj,varargin)
            %SHOW   Visualize the shape as a filled patch
            %
            % show(obj) plots the shape as a filled region using patches.
            %
            % % Examples:
            %
            % % Example1: 
            % % Create a box and visualize the structure
            % % regions
            %
            % c  = shape.Box;
            % show(c);
            % % Example 2:
            % % Create a box and plot the structure. Change the color of the
            % % box and plot it by passing a particular figure handle
            % c = shape.Box;
            % show(c);
            % c.Color = 'g';
            % fig = figure;
            % show(c,'Parent',fig,'EnableIndividualColors',true);

            % Inputs are shape(obj)
            % varargin include name value pairs:
            % EnableIndividualColors - For enabling the colors of the given
            % shape
            % Parent - Accepts figure or axes handle
            

            % Parse Inputs
            [ax,~,colflag] = parseShowInputs(obj,varargin{:});
            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end

            if getHasStructureChanged(obj)
                %           if isempty(obj.Points)
                createGeometry(obj);
            end
            updateGeometryCache(obj);
            geometry        = getGeometry(obj);
            %  antennaColor    = [223, 185, 58]/259;
            if iscell(geometry)
                BorderVertices = cell2mat(cellfun(@(x) x.BorderVertices,...
                    geometry', 'UniformOutput',false));
                [~, mult, unit] = engunits(BorderVertices);
                for m = 1:numel(obj.Geometry)
                    geometry{m}.multiplier = mult;
                    %                     geometry{m}.unit       = unit;
                end
            else
                BorderVertices = geometry.BorderVertices;
                [~, mult, unit] = engunits(BorderVertices);
                geometry.multiplier = mult;
                geometry.unit       = unit;
            end
            ScaledPoints = obj.Points*mult;

            % Unique shape id to know how many shapes are present
            usid=unique(obj.MultiShapeID);
            % BoundaryEdges
            warnState = warning('off','MATLAB:triangulation:PtsNotInTriWarnId');
            trbound=triangulation(obj.Triangles(:,1:3),ScaledPoints);
            warning(warnState);
            % trbound=em.internal.authoring.walkTriangulation(trbound);
            if isa(obj,'shape.Sphere')
                boundedges = geometry.BoundaryEdges;
                if iscell(boundedges)
                    % when scale(obj,fact,'Parent',gcf) is used
                    boundedges = cell2mat(boundedges);
                end
            else
                boundedges = featureEdges(trbound,pi/6);
            end
            %            boundedges=freeBoundary(trbound);

            % Update the patch based on the inputs given. Conditions will
            % vary depending upon how many input among
            % color,transparency,edgecolor are given and also based on
            % whether operations is performed on the object or not

            [fColor,colid] = getShowParam(obj,colflag);

            % object show before performing operation
            if isempty(obj.MultiShapeID)
                [~] = patch(ax,'Vertices',ScaledPoints,'Faces',obj.Triangles,...
                    'FaceColor',fColor,'EdgeColor','None',...
                    'faceAlpha',obj.Transparency,'EdgeAlpha',obj.Transparency);
                hold on;
                [~] = patch(ax,'Vertices',ScaledPoints,'Faces',boundedges,...
                    'EdgeColor',obj.EdgeColor);
            else
                % Object show after performing operation
                for i = 1:size(usid,1)
                    if ~isempty(obj.ColorID)
                        patch(ax,'Vertices',ScaledPoints,'Faces',...
                            obj.Triangles(find((obj.MultiShapeID(:,1))==i),1:3),...
                            'FaceColor',colid{i,1},...
                            'EdgeColor','None','faceAlpha',...
                            colid{i,2},'EdgeAlpha',colid{i,2});
                    else
                        % No boolean but multishapeid is present because of
                        % mesh (when object display used after mesh)
                        patch(ax,'Vertices',ScaledPoints,'Faces',...
                            obj.Triangles(find((obj.MultiShapeID(:,1))==i),1:3),...
                            'FaceColor','y',...
                            'EdgeColor','k','faceAlpha',...
                            1,'EdgeAlpha',1);
                    end
                    hold on;
                end
                if ~isempty(obj.ColorID)
                    [~]=patch(ax,'Vertices',ScaledPoints,'Faces',boundedges,'EdgeColor',obj.ColorID{i,3});
                end
            end
            
            if ~isempty(obj.Tetrahedra)
                [~]=tetramesh(obj.Tetrahedra,'FaceColor','y');
            end
            % Added for fixing missing Z-axis ticks g2722086
            X= ScaledPoints(:,1);
            Y= ScaledPoints(:,2);
            Z= ScaledPoints(:,3);
            axis(ax,'equal');
            em.MeshGeometry.decoratefigureandaxes(ax, X, Y, Z, unit,[57 9]);
            % Added for image going out of figure window g2722173
            box(ax,'on'); grid(ax,'on');     
        end



        function varargout = scale(obj,scaling,varargin)
            %SCALE   Change the size of the shape by a fixed amount
            %
            % OBJ = scale(obj,scaling) scales the shape by a constant factor
            % % Examples:
            % 
            % % Example 1:
            %
            % % Create a box and scale it by a factor 2
            % sh = shape.Box;
            % scale(sh,2);

            % inputs are shape (obj) ,scaling factor(scaling)
            % varargin include Reference point(RefPoint),name - value pair 
            % Parent 
           

            if isempty(obj)
                error(message("antenna:shapeerrors:EmptyShape"));
            end
            if nargin==2
                RefPoint =[];
            end
            if nargin ==3 ||nargin==5
                RefPoint = varargin{1};
                validateattributes(RefPoint,{'numeric'},...
                    {'nonempty','nonnan','finite','real',...
                    'numel',3},'scale','RefPoint',3);

            end
            %             RefPoint=varargin{1};
            if nargin==4
                RefPoint=[];
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{1:2});
                ParentObj=parserObj.Results.Parent;
                checkParentInput(obj,ParentObj);

            elseif nargin ==5
                parserObj=inputParser;
                addParameter(parserObj,'Parent',[]);
                parse(parserObj,varargin{2:3});
                ParentObj=parserObj.Results.Parent;
                checkParentInput(obj,ParentObj);

            else
                ParentObj=[];
            end


            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            if isa(obj,'shape.Circle')&& ~(isscalar(scaling))
                error(message('antenna:shapeerrors:OnlyScalar', ...
                    'scaling for shape.Circle'));
            end
            % Scale
            S = scaling;
            %             RefPoint = obj.Center;

            vert = getShapeVertices(obj);
            [xcenter]  = (max(vert(:,1))+min(vert(:,1)))/2;
            [ycenter]  = (max(vert(:,2))+min(vert(:,2)))/2;
            [zcenter]  = (max(vert(:,3))+min(vert(:,3)))/2;

            if ~isempty(RefPoint)
                % Check if the shape is symmetrical about x and y axis.
                % The scale operation, with respect to the centroid is valid for
                % all the symmetrical shapes as the centroid of the
                % shape is the center. Hence the reference point has to
                % be ignored if the shape is symmetrical and consider
                % if it is not.
                V = getShapeVertices(obj);
                V = unique(V,'rows','stable');
                A = V(:,1:3);
                % sort the y points
                A2y = unique((sort(A(:,2),'ascend')));
                A2z = unique((sort(A(:,3),'ascend')));
                % Find x points for the sorted y points
                x1 = cell(1,size(A2y,1));
                x2 = cell(1,size(A2z,1));
                x = [];
                xx=[];
                for i = 1:size(A2y,1)
                    x1{i} = find(abs(A(:,2)-A2y(i)) <= eps(max(A(:,2),A2y(i))));
                    x(end+1:end+size(x1{i},1),:)    = x1{i};
                end
                for i = 1:size(A2z,1)
                    x2{i} = find(abs(A(:,3)-A2z(i)) <= eps(max(A(:,3),A2z(i))));
                    xx(end+1:end+size(x2{i},1),:) =x2{i};
                end

                x = intersect(x,xx,'rows');
                x2 = unique(x,'row','stable');
                Ax = A(x2,1);
                Ay = sort(A(:,2),'ascend');
                Az = sort(A(:,3),'ascend');
                % Vertices after sorting
                A2 = [Ax,Ay,Az];
                % check the distances from the centroid to all the points
                distances = round(sqrt((A2(:,1)-xcenter).^2 + (A2(:,2)-ycenter).^2+(A2(:,3)-zcenter).^2),12);
                n = size(A2,2);
                if mod(n,2)==1
                    n1 = round(floor(n/2),12);
                else
                    n1 = n/2;
                end
                % Distance from the top to the mid points and from the bottom to the mid
                % points should be same if the structure is symmetry
                d1 = distances(1:n1);
                d1 = sort(d1,'ascend');
                d2 = flip(distances);
                d3 = d2(1:n1);
                d3 = sort(d3,'ascend');
                if all(d1==d3)||isprop(obj,'Center')
                    warning(message('antenna:shapeerrors:IgnoreRefPoint'))
                    center_x =  xcenter;
                    center_y =  ycenter;
                    center_z =zcenter;
                else
                    center_x =  RefPoint(1);
                    center_y =  RefPoint(2);
                    center_z =RefPoint(3);
                end
            else
                center_x =  xcenter;
                center_y =  ycenter;
                center_z = zcenter;
            end

            % Reference point
            RefPoint =[ center_x center_y center_z];
            % Shift reference to origin,perform scaling w.r.to origin,move it back
            newpt = obj.Points-RefPoint;
            newpt = newpt.*S;
            newpt =newpt+RefPoint;

            % Assign scaled points
            obj.Points=newpt;

            updateShapeDimension(obj,S);
            resetHasStructureChanged(obj);
            % 3. Update Mesh if it exists
            isMeshed = ~isempty(obj.MesherStruct.Mesh.p);

            if isMeshed
                obj.MesherStruct.Mesh.p = (S*obj.MesherStruct.Mesh.p);
                obj.Points = obj.MesherStruct.Mesh.p;
            end

            if nargout==0
                if isempty(ParentObj)
                    show(obj);
                else
                    show(obj,"Parent",ParentObj);
                end
            else
                varargout{1} = obj;
            end
        end
        
    end
    methods(Hidden)
        function rotate3D(obj,ang,axispt1,axispt2)

            % Rotates the vertices by an angle(ang) around axis defined by
            % axispt1,axispt2
            
            % 1. Update BorderVertices
            rotateVertices(obj,ang,axispt1,axispt2);
            % Get the tilt value and assign it to tilt
            resetState = getHasStructureChanged(obj);
            obj.Tilt = [obj.Tilt ang];
            tmpAxis = [axispt1;axispt2];
            obj.TiltAxis = [obj.TiltAxis; tmpAxis];
            obj.MesherStruct.HasStructureChanged = resetState;

            if isprop(obj,'Center')
                rotCenter = em.internal.rotateshape(obj.Center',axispt1,axispt2,ang);
                if ~isequal(obj.Center,rotCenter')
                    obj.Center = rotCenter';
                    resetHasStructureChanged(obj);
                end
            end

            % 2. Update Mesh if it exists
            isMeshed = ~isempty(obj.MesherStruct.Mesh.p);
            if isMeshed  %--> allow rotation of mesh, show
                obj.MesherStruct.Mesh.p = em.internal.rotateshape(obj.MesherStruct.Mesh.p,axispt1,axispt2,ang);
                obj.Points = obj.MesherStruct.Mesh.p;
                obj.Points=obj.Points';
            end
        end
        function rotateVertices(obj,ang,axispt1,axispt2)
            % Rotates the vertices by an angle

            Geometry                = getGeometry(obj);
            % If geometry has no Border vertices,create geometry and assign it
            if isempty(Geometry.BorderVertices)
                createGeometry(obj);
                Geometry                = getGeometry(obj);
            end
            % Assign it to geometry after rotation
            obj.MesherStruct.Geometry.BorderVertices = em.internal.rotateshape(Geometry.BorderVertices'...
                ,axispt1,axispt2,ang)';
            obj.Points=obj.MesherStruct.Geometry.BorderVertices;
        end
   
        function        P=getpatch(usid,inp1,inp2,inp3,inp4,inp5)
            % Get the patch from inputs

            P(:,1:usid)=0;
            for i=1:size(usid,1)
                P1=patch('Vertices',inp1,'Faces',inp2,'FaceColor',inp3,'EdgeColor',inp4,'faceAlpha',inp5,'EdgeAlpha',inp6);
                P(:,i)=P1;
                hold on;
            end
        end 
      
        function translateOutOfXYPlane(obj,vec)
            % Translate the vertices

            % 1. Update BorderVertices
            translateVertices(obj,vec)
            %             resetHasStructureChanged(obj);

            % 2. Update Mesh if it exists
            isMeshed = ~isempty(obj.MesherStruct.Mesh.p);

            if isMeshed  %--> allow rotation of mesh, show
                obj.MesherStruct.Mesh.p = em.internal.translateshape(obj.MesherStruct.Mesh.p,vec);
                obj.Points = obj.MesherStruct.Mesh.p;
                obj.Points=obj.Points';
            end
        end
        
        function translateVertices(obj,vec)
            % Translate the vertices 
            
            Geometry                = getGeometry(obj);
            if isempty(Geometry.BorderVertices)
                createGeometry(obj);
                Geometry                = getGeometry(obj);
            end
            obj.MesherStruct.Geometry.BorderVertices = em.internal.translateshape(Geometry.BorderVertices'...
                ,vec)';
            obj.Points=obj.MesherStruct.Geometry.BorderVertices;
        end
  
        function   [pp,tt]  = remapdomain(~,pIn,tIn,pOut,tOut,feedsAndVias)
            % Remap the domain numbers
            % 1. Find sample point out of each domain number in the output set
            % 2. Check where that sample point lands in the input set
            % 3. Find all the tris with that domain number and swap out the domain
            % number.
            domainNumIn     = unique(tIn(:,4));
            domainNumOut    = unique(tOut(:,4));
            if isequal(numel(domainNumIn),1)
                tOut(:,4) = domainNumIn;
            elseif ~isequal(domainNumIn,domainNumOut)
                warnflag = warning('Off','MATLAB:triangulation:PtsNotInTriWarnId');
                TRin    = triangulation(tIn(:,1:3),pIn);
                TRout   = triangulation(tOut(:,1:3),pOut);
                eIn     = edges(TRin);
                eOut    = edges(TRout);
                warning(warnflag);
                inCenters = incenter(TRin); %em.internal.meshprinting.meshtricenter(pIn, tIn);
                inNormals = faceNormal(TRin);
                outCenters = incenter(TRout); %em.internal.meshprinting.meshtricenter(pOut,tOut);
                outNormals = faceNormal(TRout);
                newDomain = tOut(:,4);
                for i = 1:numel(domainNumOut)
                    testDomain = domainNumOut(i);
                    testTri = find(tOut(:,4) == testDomain);
                    tempNormals = outNormals(testTri,:);
                    tempCenters = outCenters(testTri,:);
                    % Find distance between test triangles in out mesh and centers in
                    % input mesh and index of closest center
                    [k,dist] = dsearchn(inCenters,tempCenters);
                    [~,ind]=sort(dist);
                    k=k(ind);
                    %         k = unique(k);
                    isDomainSet = false;
                    ctr = 0;
                    % Start with closest point in k, and run through all normals in
                    % output mesh for this domain. If any are close, assign and break
                    % and move to next.
                    while ~isDomainSet && ctr < numel(k)
                        ctr = ctr + 1;
                        % Checking normals
                        for q = 1:size(tempNormals,1)
                            if ~isempty(k) && abs(abs(dot(tempNormals(q,:),inNormals(k(ctr),:))))>0.5
                                newDomain(testTri) = tIn(k(ctr),4);
                                isDomainSet = true;
                                break;
                            end
                        end
                    end

                end
                tOut(:,4) = newDomain;
                % Re-check input and output domains
                if ~isequal(domainNumIn,unique(newDomain))
                    % Check the feed edge triangles
                    if ~isempty(feedsAndVias)
                        for i = 1:size(feedsAndVias,1)
                            [Dout, INDout] = em.internal.meshprinting.inter2_point_seg(pOut,eOut,feedsAndVias(i,:) );
                            [Din, INDin] = em.internal.meshprinting.inter2_point_seg(pIn,eIn,feedsAndVias(i,:) );
                            tempIndex = find(Dout<sqrt(eps) & INDout==-1);
                            if ~isempty(tempIndex)&& isscalar(tempIndex)
                                indexOut(i) = tempIndex; %#ok<AGROW>
                            else
                                closestEdgeId = find(Dout==min(Dout));
                                indexOut(i) = closestEdgeId(1); %#ok<AGROW>
                            end
                            tempIndex = find(Din<sqrt(eps) & INDin==-1); 
                            if ~isempty(tempIndex)&& isscalar(tempIndex)
                                indexIn(i) = tempIndex; %#ok<AGROW>
                            else
                                closestEdgeId = find(Din==min(Din));
                                indexIn(i) = closestEdgeId(1); %#ok<AGROW>
                            end
                            % Fixing the domains for the triangles attached to the feed
                            % edge
                            feed_edges = eOut(indexOut,:);
                            ID_outTri = edgeAttachments(TRout,feed_edges);
                            ID_outTri = cell2mat(ID_outTri);
                            feed_edges = eIn(indexIn,:);
                            ID_inTri = edgeAttachments(TRin,feed_edges);
                            ID_inTri = cell2mat(ID_inTri);
                            tOut(ID_outTri,4) = tIn(ID_inTri,4);
                        end
                    end
                end
            end

            % Check if multiple input domains have any 0's and increase by 1
            if numel(domainNumIn) > 1
                if any(domainNumIn==0)
                    tOut(:,4) = tOut(:,4) + 1;
                end
            end


            pp = pOut';
            tt = tOut';

        end
     
    end
    methods(Access=protected)

        function updateShapeDimension(obj, ~)
            if isa(obj,'shape.Custom3D')
                return;
            end
        end
    end
    methods
        function  vert = getShapeVertices(obj)
            %             if isempty(obj.Points)
            if getHasStructureChanged(obj)
                createGeometry(obj);
            end
            vert = obj.Points;
        end
    end
    methods(Hidden)
        function rObj = superLoadShape(obj,s)
            % if em.internal.authoring.custom3DFeature('EnableAuthoring')
                obj.Triangles =  s.Triangles;
                obj.Points = s.Points;
                obj.MultiShapeID =s.MultiShapeID;
                obj.Metal=s.Metal;
                obj.Dielectric = s.Dielectric;
                obj.Color=s.Color;
                obj.Transparency=s.Transparency;
                obj.EdgeColor =s.EdgeColor;
                obj.MaterialID = s.MaterialID;
                obj.TetID   = s.TetID;
                obj.ColorID = s.ColorID;
                obj.HasMaterialID = s.HasMaterialID;
                obj.HasTetID   = s.HasTetID;
                rObj = obj;
            % else
            %     error(message('antenna:shapeerrors:UnsupportedFunctionality'));
            % end
        end
        function s = saveobj(obj)           %---SAVEOBJ IMPLEMENTATION----%
            obj.MesherStruct.CacheFlag = obj.MesherStruct.HasStructureChanged;
            s = obj;
        end
        function onlyDiel = isOnlyDielectric(shape)
            onlyDiel = (isempty(shape.TetID)&&~strcmpi(shape.Dielectric,'Air'))...
                      ||(~isempty(shape.TetID)&&all(cellfun(@(x) ~strcmpi(x,'Air'),shape.TetID(:,1))));
        end
    end
    methods(Static,Hidden)
        function r = loadobj(obj)
            % checkFeatureFlag(obj);
            r = obj;
            r.MesherStruct.HasStructureChanged = obj.MesherStruct.CacheFlag;
        end
    end
end
