classdef gear < handle
    %GEAR Create gear object
    %
    % Creation:
    %   g = gear(gearRack(m,__),z)
    %   g = gear(__,'-x',x)
    %   g = gear(__,'-beta',beta)
    %   g = gera(__.'-np',np)
    %
    % Methods:
    %   print(g,[fid])
    %   plot(g,varargin)
    %       plot(__,'-fig',figno)
    %       plot(__,{'keep','nocir','nocen','notxt','save'.'gen'})
    %           'keep' keep old figure
    %           'nocir' do not plot gear characteristic circles
    %           'nocen' do not mark gear ceneter point
    %           'notxt' do not plot title
    %           'save'  save figure as jpg file
    %           'gen'   plot gear involute formation   
    %   [X,Y] = gearContour(obj,'-file',fnam) % export gear contour as XYZ
    %   file that can be imported to SOLDIWORKS
    %
    % Involute data: t=0..1 (normalized coordinate 0=top, 1=bottom)
    %   t = getPar(obj,R)  ... obtain normalized parameter for given radius
    %   s = getData(obj,t)
    %       s.X,s.Y - point on involute
    %       s.R     - radius
    %       s.L     - =sqrt(R^2 - Rb^2) - roll length
    %       s.phi   - angle in degrees = atan(X(Y)
    %       s.s     - tooth thickness
    %       s.th    - involute angle = acos(Rg/R)
    %       s.invth - involute of th in degrees = tan(th) - th
    %   [x,y,nx,ny] = getNormal(obj,t)
    %   [x,y,tx,ty] = getTangent(obj,t)
    %   [x,y,Fx,Fy,M] = getForce(obj,t)
    %
    % Other methods:
    %   plotRackContact(g)
    
    properties (Dependent)
        z          % number of teeth
        x          % profile shift coefficient
        beta       % helix angle in degrees
        u          % Tip shortening coefficient
    end
    properties (Access = private)
        z_          % number of teeth
        x_          % profile shift coefficient
        beta_       % helix angle in degrees
        u_          % Tip shortening coefficient
    end
    
    % not used
    properties (Dependent)
        Ri          % gear inner radius
    end
    properties (Access = private)
        Ri_         % gear inner radius
    end
    
    properties (Dependent)
        np   % number of points on tooth profile
    end
    properties (Access = private)
        np_
    end
    
    properties
        rack        % basic rack, i.e., gearRack object
    end
    
    properties(Dependent)
        R0   % reference pitch circle radius
        Rc   % max. radius for pointed teeth        
        Ra   % addendum (tip) circle radius
        Rr   % pitch circle radius
        Ru   % undercut circle
        Rb   % base circle radius
        Rd   % dedendum (root) circle radius
        e    % profile shift
        h    % total tooth height
        ha   % addendum tooth height
        hd   % dedendum tooth height
        sa   % tooth thickness at the addendum circle
        sr   % tooth thickness at the reference pitch circle
        su   % tooth thickness at undercuting circle
        sb   % tooth thickness at the base circle
        sd   % tooth thicknes at root circle
        xmax % max. profile shift coefficient
        xmin % min. profile shift coefficient
        zmin % minimum nuber of teeth to avoud undercutting
        L    % length of segments
        Lt   % total length of profile
    end
    
    properties (Dependent)
        tpar % end point values of parameter
        tc   % max. radius parameter
    end
    
    properties (Access = private)
        ts_   % start parameter value
        te_   % end parameter value
        Rc_   % intersection with y-axis (ultimate gear outer radius)
        xmax_ % max shift to achive pointed teeth
        tc_
    end
    
    methods
        function obj = gear(rack,z,varargin)
            
            narginchk(2,inf)
            
            % check input
            validateattributes(z,  {'numeric'}, {'>',1,'integer','scalar'});
            if ~isa(rack,'gearRack')
                error('expect gearRack data.')
            end
            
            % set default values
            beta = 0;
            x    = 0;
            np   = 100;
            
            % scann optional input name-value pairs
            if ~isempty(varargin)
                for k = 1:2:length(varargin)
                    switch lower(varargin{k})
                        case '-beta'
                            beta = varargin{k + 1};
                            validateattributes(beta,  {'numeric'}, {'>=',0,'<',90,'real','scalar'});
                        case '-x'
                            x = varargin{k + 1};
                            validateattributes(x,    {'numeric'}, {'real','scalar'});
                        case '-np'
                            np = varargin{k + 1};
                            validateattributes(np,     {'numeric'}, {'>',3,'<',101,'integer','scalar'});
                        otherwise
                    end
                end
            end
            
            % save data
            obj.rack = rack;
            obj.z_ = z;
            obj.x_ = x;
            obj.beta_ = beta;
            obj.np_ = np;
            obj.u_ = obj.rack.u;
            obj.Ri_ = 0;
            tinit(obj);
        end
    end
    
    % Set methods
    methods
        function set.beta(obj,beta)
            % set helix angle in degree
            validateattributes(beta,{'numeric'},{'>=',0,'<',90,'real','scalar'});
            obj.beta_ = beta;
            tinit(obj)
        end
        function set.np(obj,np)
            % set number of points on the profile
            validateattributes(np,{'numeric'}, {'>',9,'<',201,'integer','scalar'});
            obj.np_ = np;
        end
        function set.Ri(obj,Ri)
            % set inner radius
            validateattributes(Ri,{'numeric'}, {'>=',0,'<',obj.Rd,'real','scalar'});
            obj.Ri_ = Ri;
        end
        function set.x(obj,x)
            % set profile shift coefficient
            validateattributes(x,{'numeric'}, {'real','scalar'});
            obj.x_ = x;
            tinit(obj)
        end
        function set.z(obj,z)
            % set number of teeth
            validateattributes(z,  {'numeric'}, {'>',1,'integer','scalar'});
            obj.z_ = z;
            tinit(obj)
        end
        function set.u(obj,u)
            % set tip shortening coefficient
            validateattributes(u,  {'numeric'}, {'>',0,'<=',1,'real','scalar'});
            obj.u_ = u;
            obj.rack.u = u;  % change rack data !!!
            tinit(obj)
        end
    end
    
    methods
        function tinit(obj)
            % calculate key points parameters only when data are changed
            [ts,te] = calcParam(obj);
            if te(1) < ts(1)
                te(1) = ts(1);  % avoid pointed teeth
            end
            obj.ts_ = ts;
            obj.te_ = te;
            [obj.Rc_,obj.tc_] = calcRc(obj);
            
            if obj.tc_ > ts(2)
                obj.ts_(2) = obj.tc_;
            end
            obj.xmax_ = calcXmax(obj);
        end
    end
    
    % Get methods for data
    methods
        function val = get.beta(obj)
            % helix angle
            val = obj.beta_;
        end
        function val = get.Ri(obj)
            % inner radius
            val = obj.Ri_;
        end
        function val = get.x(obj)
            % profile shift factor
            val = obj.x_;
        end
        function val = get.z(obj)
            % number of teeth
            val = obj.z_;
        end
        function val = get.u(obj)
            % tooth shortening coefficient
            val = obj.u_;
        end
        function val = get.np(obj)
            % number of points for contour
            val = obj.np_;
        end
    end
    
    % Get methods for dependent properties
    methods
        function val = get.e(obj)
            % profile shift
            val = obj.x*obj.rack.m;
        end
        function val = get.h(obj)
            % total tooth height
            val = obj.Ra - obj.Rd;
        end
        function val = get.ha(obj)
            % addendum tooth height
            if obj.te_(1) == 0
                val = obj.Rc_ - obj.Rr;
                return
            end
            val = obj.Ra - obj.Rr;
        end
        function val = get.hd(obj)
            % dedendum tooth height
            val = obj.Rr - obj.Ru;
        end
        function val = get.L(obj)
            % length of profile segments
            val = calcLength(obj);
        end
        function val = get.Lt(obj)
            % total length  of profile
            val = sum(obj.L);
        end
        function val = get.Rc(obj)
            % max. circle
            val = obj.Rc_;
        end
        function val = get.Ru(obj)
            % undercut circle
            %  [ts,~] = calcParam(obj);
            [X,Y] = profil(obj,3,obj.ts_(3));
            val = sqrt(X^2 + Y^2);
        end
        function val = get.R0(obj)
            % reference pitch circle
            val = obj.z*obj.rack.m/2;
        end
        function val = get.Rr(obj)
            %  pitch circle
            val = obj.R0/cosd(obj.beta);
        end
        function val = get.Ra(obj)
            % addendum (tip) circle
            val = obj.Rr + obj.rack.m*(obj.rack.u + obj.x);
            if val > obj.Rc_
               % val = obj.Rc_;
            end
        end
        function val = get.Rb(obj)
            % base circle
            val = obj.R0*cosd(obj.rack.alpha)/sqrt(1 -(sind(obj.beta)*cosd(obj.rack.alpha))^2);
        end
        function val = get.Rd(obj)
            % dedendum (root) circle
            val = obj.R0/cosd(obj.beta) - obj.rack.m*(1 + obj.rack.c - obj.x);
        end
        function val = get.sa(obj)
            % tooth thickness at root circle
            if obj.te_(1) == 0
                val = 0;
                return
            end
            val = toothThickness(obj,obj.Ra);
        end
        function val = get.sr(obj)
            % tooth thickness at pitch circle
            val = toothThickness(obj,obj.Rr);
        end
        function val = get.su(obj)
            % tooth thickness at involute start circle
            val = toothThickness(obj,obj.Ru);
        end
        function val = get.sb(obj)
            % tooth thickness at base circle
            val = toothThickness(obj,obj.Rb);
        end
        function val = get.sd(obj)
            % tooth thickness at root circle
            %val = toothThickness1(obj,obj.Rd);
            %return
            [~,~,rho,~,~] = calcKeyPoints(obj.rack);
            val = 2*obj.Rd/obj.R0*(obj.rack.m*(pi/4 + tand(obj.rack.alpha)) + rho*cosd(obj.rack.alpha));
        end
        function val = get.xmax(obj)
            % max. profile shif coefficient
            val = obj.xmax_;            
        end
        function val = get.xmin(obj)
            % min profile shift coefficient to avoud undercutting
            val = 1 - obj.z/obj.zmin;
        end
        function val = get.zmin(obj)
            % min. number of teeth to avoid undercuting
            sina = sind(obj.rack.alpha);
            cosb = cosd(obj.beta);
            cosa = sqrt(1 - sina^2);
            sinb = sqrt(1 - cosb^2);
            val = 2*cosb*(1 - (sinb*cosa)^2)/sina^2;
        end
        function out = get.tpar(obj)
            % parameters for key points
            out = [obj.ts_' obj.te_'];
        end
        function out = get.tc(obj)
            % parameters for key points
            out = obj.tc_;
        end        
    end
    
    % involute segment get function
    methods
        function t = getPar(obj,R)
            % calculate normalized parameter for given radius R (involute segment
            % only)
            validateattributes(R,{'numeric'},{'>=',obj.Ru,'<=',obj.Ra,'real','scalar'});           
            try
                %==============================
                t = mfzero(@fun,[obj.ts_(2),obj.te_(2)]);
                %===============================
                t = (t - obj.ts_(2))/(obj.te_(2) - obj.ts_(2));
            catch
                t = NaN;
            end            
            function val = fun(t)
                [X,Y] =  profil(obj,2,t);
                val = X^2 + Y^2 - R^2;
            end
        end        
        function s = getData(obj,t)
            % get various geometric data of involute at point given
            % by parameter t
            narginchk(2,2)
            % check input
            validateattributes(t,{'numeric'}, {'>=',0,'<=',1,'real','scalar'});
            t = obj.ts_(2) + t*(obj.te_(2) - obj.ts_(2));
            [X,Y] = profil(obj,2,t);
            s.t   = t;
            s.R   = sqrt(X^2 + Y^2);  % radius            
            s.X   = X;
            s.Y   = Y;
            s.phi = atand(X/Y);  % angle in degree
            s.s   = s.R*(2*s.phi)*pi/180; % thickness
            s.th  = acosd(obj.Rb/s.R);  % involute angle
            s.invth = tand(s.th)*180/pi - s.th;  % involute of th
            s.L   = sqrt(s.R^2 - obj.Rb^2);  % distanc between tangent point point on involute            
        end
        function [X,Y,Tx,Ty] = getTangent(obj,t)
            % get point and tangent to involute point given with parameter t
             narginchk(2,2)
             nargoutchk(4,4)            
            % check input
            validateattributes(t,{'numeric'}, {'>=',0,'<=',1,'real','scalar'});
            t = obj.ts_(2) + t*(obj.te_(2) - obj.ts_(2));
            [X,Y] = profil(obj,2,t);
            [Tx,Ty] = calcTangent(obj,2,t);
        end        
        function [X,Y,Nx,Ny] = getNormal(obj,t)
            % get point and normal to involute point given with parameter t
             narginchk(2,2)
             nargoutchk(4,4)            
            % check input
            validateattributes(t,{'numeric'}, {'>=',0,'<=',1,'real','scalar'});
            t = obj.ts_(2) + t*(obj.te_(2) - obj.ts_(2));
            [X,Y] = profil(obj,2,t);
            [dX,dY] = calcTangent(obj,2,t);
            Nx = -dY;
            Ny = dX;
        end
        function [X,Y,Fx,Fy,M] = getForce(obj,t,F)
            % get components of force at given point on involute and its
            % moment with respect to center point
            [X,Y,Nx,Ny] = getNormal(obj,t);
            Fx = F*Nx;
            Fy = F*Ny;
            M  = X*Fy - Y*Fx; 
        end
    end
    
    % Output methods
    methods
        
        function plot(obj,varargin)
            % Plot gear
            % set default values
            fig = [];
            nz = obj.z;
            init = true;
            drwcirc = true;
            drwcent = true;
            drwtxt = true;
            save = false;
            gen = false; % generation of tooth
            % scann optional input name-value pairs
            if ~isempty(varargin)
                id = zeros(length(varargin),1); % use for delition of options
                for k = 1:length(varargin)
                    if ~ischar(varargin{k})
                        continue
                    end
                    switch lower(varargin{k})
                        case '-fig'
                            fig = varargin{k + 1};
                            validateattributes(fig,  {'numeric'}, {'positive','integer','scalar'});
                            id(k:k+1) = 1;
                        case {'gen','gentooth','generationoftooth'}
                            gen = true;
                            id(k) = 1;
                        case '-nz'
                            nz = varargin{k + 1};
                            validateattributes(nz,{'numeric'}, {'>',0,'<=',obj.z,'integer','scalar'});
                            id(k:k+1) = 1;
                        case 'keep'
                            init = false;
                            id(k) = 1;
                        case {'nocir','nocirc'}
                            drwcirc = false;
                            id(k) = 1;
                        case {'nocent','nocenter','nocen'}
                            drwcent = false;
                            id(k) = 1;
                        case {'notxt','notext'}
                            drwtxt = false;
                            id(k) = 1;
                        case 'save'
                            save = true;
                            id(k) = 1;
                        otherwise
                    end
                end
                % delete used options
                varargin(id == 1) = [];
            end
            % plot begin
            if init
                if isempty(fig)
                    drawInit
                else
                    drawInit(fig)
                end
                axis off
            else
                hold on
            end
            
            if gen
                plotEnvelope(obj)
            else
                if drwcirc
                    % add characteristic circles
                    drawCircle(0,0,obj.Rr,'r:')
                    drawCircle(0,0,obj.Rb,'g:')
                    drawCircle(0,0,obj.Rd,'b:')
                    drawCircle(0,0,obj.Ra,'m:')
                    drawCircle(0,0,obj.Ru,'k:')
                end
                if drwcent
                    %plot center point
                    drawPoint(1,obj.rack.m/4,0,0)
                end
                if drwtxt
                    % add title
                    txt = sprintf('m = %g  z = %g  %s = %g%s  d = %g  d%s = %g  x = %g  h = %g',...
                        obj.rack.m,obj.z,'\beta',obj.beta,'^0',2*obj.R0,'_b',2*obj.Rb,obj.x,obj.h);
                    title(txt,'FontSize',12,'FontWeight','normal')
                end
                % plot gear
                [X,Y] = calcPoints(obj,nz);
                if ~isempty(varargin)
                    drawPolyline(X,Y,varargin{:})
                else
                    drawPolyline(X,Y,'k','LineWidth',2)
                end
                % draw inner circle
                if obj.Ri > 0
                    if ~isempty(varargin)
                        drawCircle(0,0,obj.Ri,varargin{:})
                    else
                        drawCircle(0,0,obj.Ri,'k','LineWidth',2)
                    end
                end
                % limits for plot
                if init
                    drawLimits(min(X),max(X),min(Y),1.1*max(Y));
                end
            end
            if save
                fnam = drawSave;
                fprintf('Drawing is saved to the file %s\n',fnam);
            end
        end
        
        function  print(obj,fid)
            if nargin < 2
                fid = 1;
            end
            [ad,am,as]=deg2dms(obj.rack.alpha);
            [bd,bm,bs]=deg2dms(obj.beta);
            fprintf(fid,'Data\n');
            fprintf(fid,'                               Number of teeth:%12d\n',obj.z);
            fprintf(fid,'                                        Module:%12.4f mm\n',obj.rack.m);
            fprintf(fid,'                       Standard pressure angle:%6d:%02d:%02d \n',ad,am,as);
            fprintf(fid,'             Manufacturing tip tooth clearance:%12.4f mm\n',obj.rack.c);
            fprintf(fid,'                The tip shortening coefficient:%12.4f \n',obj.u);
            fprintf(fid,'                                   Helix angle:%6d:%02d:%02d \n',bd,bm,bs);
            fprintf(fid,'                    Profile shift coefficients:%12.4f\n',obj.x);
            fprintf(fid,'Calculated parameters\n');
            fprintf(fid,'                Addendum (tip) circle diameter:%12.4f mm\n',2*obj.Ra);
            fprintf(fid,'               Reference pitch circle diameter:%12.4f mm\n',2*obj.Rr);
            fprintf(fid,'                      Undercut circle diameter:%12.4f mm\n',2*obj.Ru);
            fprintf(fid,'                          Base circle diameter:%12.4f mm\n',2*obj.Rb);
            fprintf(fid,'               Dedendum (root) circle diameter:%12.4f mm\n',2*obj.Rd);
            fprintf(fid,'                            Total tooth height:%12.4f mm\n',obj.h);
            fprintf(fid,'                         Addendum tooth height:%12.4f mm\n',obj.ha);
            fprintf(fid,'                         Dedendum tooth height:%12.4f mm\n',obj.hd);
            fprintf(fid,'        Tooth thickness at the addendum circle:%12.4f mm\n',obj.sa);
            fprintf(fid,'           Tooth thickness at the pitch circle:%12.4f mm\n',obj.sr);
            fprintf(fid,'        Tooth thickness at the undercut circle:%12.4f mm\n',obj.su);
            fprintf(fid,'            Tooth thickness at the base circle:%12.4f mm\n',obj.sb);
            fprintf(fid,'            Tooth thickness at the root circle:%12.4f mm\n',obj.sd);
            if obj.Ru > obj.Rb && obj.z < obj.zmin
                fprintf(fid,'***Warning: gear is undercut.\n');
            end
            if obj.sa < 0.25*obj.sr
                fprintf(fid,'***Warning: pointed teeth.\n');
            end
        end
    end
    
    methods
        %=======================
        % Aux. plot functions
        %========================
        function plotNormal(obj,t,t1,t2,varargin)
            % plot normal at given point on involute 
            [X,Y,Nx,Ny] = getNormal(obj,t);
            for k = 1:length(t)
                drawLine(X(k),Y(k),X(k) + Nx,Y(k) + Ny,t1,t2,varargin{:})
            end
        end
        function plotTangent(obj,t,t1,t2,varargin)
            % plot normal at given point on involute
            [X,Y,Tx,Ty] = getTangent(obj,t);
            for k = 1:length(t)
                drawLine(X(k),Y(k),X(k) + Tx,Y(k) + Ty,t1,t2,varargin{:})
            end
        end
        function plotForce(obj,t,F,varargin)
            % plot force at given point on involute
            [X,Y,Nx,Ny] = getNormal(obj,t);
            ad1 = obj.h/6;
            ad2 = ad1/2;
            %plot from out to in
            drawArrow(-3,ad1,ad2,X,Y,'-delta',-F*Nx,-F*Ny,varargin{:})
        end
        function plotRackContactPts(obj,fig)
            % plot points on the rack that forms the tooth profile
            if nargin < 2
                drawInit
            else
                drawInit(fig)
            end
            axis off
            title('Contact sections','FontSize',12,'FontWeight','normal')
            % [ts,te] = calcParam(obj);
            [X,Y] = calcPoints(obj.rack,1,0,0,obj.np);
            drawPolyline(X,Y,'k','LineWidth',1)
            for n = 1:4
                if n ~= 3
                    [X,Y] = calcProfile(obj.rack,n,[obj.ts_(n),obj.te_(n)]);
                else
                    [X,Y] = calcProfile(obj.rack,n,linspace(obj.ts_(n),obj.te_(n),obj.np));
                end
                drawPolyline( X,Y,'k','LineWidth',4)
                drawPolyline(-X,Y,'k','LineWidth',4)
            end
        end
    end
        
    methods
        function [X,Y] = gearContour(obj,varargin)
            % ifv file is given then export contour to txt format. 
            %   This can be inported to SOLIDWORKS
            % set default values
            nz = obj.z;
            fnam = '';
            ds = obj.Lt/obj.np_;  % size step
            %dsmin = min(obj.Lt);
            %if dsmin <= 0
            %    dsmin = min(obj.Lt(2:4));
            %end
            %ds = dsmin/2
            % scan options
            if ~isempty(varargin)
                for k = 1:2:length(varargin)
                    switch lower(varargin{k})
                        case '-ds'
                            ds = varargin{k + 1};
                            validateattributes(ds,{'numeric'}, {'>',0,'<=',obj.Lt/4,'real','scalar'});
                        case {'-fnam','-fname','-file','-filename'}
                            fnam = varargin{k + 1};
                            validateattributes(fnam,{'char'},{'nonempty'});
                        case {'-np','-nn'}                            
                            nn = varargin{k + 1};
                            validateattributes(nn,{'numeric'}, {'>',3,'<=',100,'integer','scalar'});
                            ds = obj.Lt/nn;
                    end
                end
            end
            
            % Boundary points
            [X,Y] = calcPoints1(obj,ds,nz);
            if isrow(X)
                X = X';
            end
            if isrow(Y)
                Y = Y';
            end
            % close contour
            X(end+1) = X(1);
            Y(end+1) = Y(1);
            
            if isempty(fnam)
                return
            end
            
            fid = fopen(fnam,'w');
            %  fprintf(fid,'%12s%12s%12s\n','x','y','z')
            for n = 1:length(X)
                %fprintf(fid,'%12.4f\t%12.4f\t%12.4f\n',node(n,1),node(n,2),0);
                fprintf(fid,'%.4g\t%.4g\t%.4g\r\n',X(n),Y(n),0);
            end
            fclose(fid);
        end
    end
    
    % FEM model
    methods
       function [node,edge,node1] = model(obj,varargin)
            %MODEL generate point model of conture and generate internal
            %points
            
            % set default values
            ds = obj.Lt/obj.np;
            nz = obj.z;
            
            % scan options
            if ~isempty(varargin)
                for k = 1:2:length(varargin)
                    switch lower(varargin{k})
                        case '-nz'
                            nz = varargin{k + 1};
                            validateattributes(nz,{'numeric'},{'>',0,'<=',obj.z,'integer','scalar'});
                        case '-ds'
                            ds = varargin{k + 1};
                            validateattributes(ds,{'numeric'}, {'>',0,'<=',obj.Lt/4,'real','scalar'});
                        otherwise
                    end
                end
            end
            
            % Boundary points
            [X,Y] = calcPoints1(obj,ds,nz);
            if isrow(X)
                X = X';
            end
            if isrow(Y)
                Y = Y';
            end
            node = [X Y];
            nnd = length(node);
            if nz < obj.z
                [X,Y] = evalLine(node(nnd,1),node(nnd,2),0,0,linspace(0,1,obj.Rd/ds)');
                X = [node(:,1);X];
                Y = [node(:,2);Y];
                [X,Y] = deleteDuplicate(X,Y);
                node = [X Y];
                [X,Y] = evalLine(0,0,node(1,1),node(1,2),linspace(0,1,obj.Rd/ds)');
                X = [node(:,1);X];
                Y = [node(:,2);Y];
                [X,Y] = deleteDuplicate(X,Y);
                node = [X Y];
            end
            nnd = length(node);
            edge = [(1:(nnd-1))' (2:nnd)'; nnd 1];
            
            % generate internal points
            if nz == obj.z
                k = 1;
                node1(1,1) = 0;  % center
                node1(1,2) = 0;
                nr = fix(obj.Rd/ds);
                dr = obj.Rd/nr;
                for r = dr:dr:obj.Rd - dr
                    s = 2*pi*r;
                    n = fix(s/ds);
                    if n < 4
                        continue
                    end
                    dth = 360/n;
                    for th = 0:dth:360-dth
                        k = k + 1;
                        node1(k,1) = r*cosd(th);
                        node1(k,2) = r*sind(th);
                    end
                end
                % tooth
                X = [];
                Y = [];
                j = 0;
                nr = fix((obj.Ra - obj.Rd)/ds);
                dr = (obj.Ra - obj.Rd)/nr;
                for r = obj.Rd :dr:obj.Ra - dr
                    s = toothThickness(obj,r)/2;
                    n = fix(s/ds);
                    if n < 1
                        continue
                    end
                    the = s/r*180/pi;
                    dth = the/n;
                    for th = -the+dth:dth:the-dth
                        j = j + 1;
                        X(j) = r*sind(th);
                        Y(j) = r*cosd(th);                        
                    end
                end
                X = X';
                Y = Y';
                XX = X;
                YY = Y;
                for n = 1:nz - 1
                    [X,Y] = trRot2d(X,Y,0,0,360/obj.z);
                    XX = [X;XX];
                    YY = [Y;YY];
                end
                X = XX;
                Y = YY;
                if isrow(X)
                    X = X';
                end
                if isrow(Y)
                    Y = Y';
                end
                [X,Y] = deleteDuplicate(X,Y);
                node2 = [X Y];
                node1 = [node1; node2];
            else
               % fprintf('Ups\n')
            end
            %{
             return   
             % generate internal points
             ths = atan2d(node(1,2),node(1,1));
             the = atan2d(node(nnd,2),node(nnd,1));
             if the < ths
                 the = the + 360;
             end
             dr = ds;
             X = [];
             Y = [];
             for r = dr:dr:obj.Rd-dr
                 s = r*(the - ths)*pi/180;
                 npt = fix(s/ds) + 1;
                 dth = (the - ths)/npt;
                 for j = 2:npt
                     k = k + 1;
                     X(k) = r*cosd(j*dth);
                     Y(k) = r*sind(j*dth);
                 end
             end
             if isrow(X)
                 X = X';
             end
             if isrow(Y)
                 Y = Y';
             end
             node1 = [X Y];
             % scatter(X,Y)
             
             else
                 node(nnd + 1,1:2) = [0 0];
        end
            % remove duplicate points
            %[~, I, ~] = unique(node,'first','rows');
            %I = sort(I);
            %node = node(I,:);
            nnd = length(node);
            edge = [(1:(nnd-1))' (2:nnd)'; nnd 1];
            %}
        end
    end
    
    methods (Access = private)
        function plotEnvelope(obj)
            % Plot formation of the tooth
            title('Generation of gear tooth','FontWeight','normal','FontSize',12)
            
            % draw characteristic circles
            phi = -180/obj.z;
            drawCircle(0,0,obj.Rr,phi+90,360/obj.z,'r:')
            drawCircle(0,0,obj.Rb,phi + 90,-2*phi,'g:')
            drawCircle(0,0,obj.Rd,phi + 90,-2*phi,'b:')
            drawCircle(0,0,obj.Ra,phi + 90,-2*phi,'m:')
            drawCircle(0,0,obj.Ru,phi + 90,-2*phi,'k:')
            
            % draw profile and label points
            c = ['r','b','g','m'];
            for n = 1:4
                [X,Y] = profil(obj,n,linspace(0,1));
                scatter(X(1),Y(1),30,'k','filled')
                drawText(X(1),Y(1),num2str(n))
                drawPolyline( X,Y,'LineWidth',2,'Color',c(n))
                drawPolyline(-X,Y,'LineWidth',2,'Color',c(n))
            end
            scatter(X(end),Y(end),30,'k','filled')
            drawText(X(end),Y(end),num2str(n+1))
            
            % draw  path of point 2
            [X,Y] = point2(obj,linspace(0,1,obj.np));
            drawPolyline( X,Y,'k','LineWidth',1)
            drawPolyline(-X,Y,'k','LineWidth',1)
            scatter(X(1),Y(1),30,'k')
            drawText(X(1),Y(1),sprintf('2'''))
            
            % draw rack
            [X,Y] = calcPoints(obj.rack, 1, 0,0);
            X0 = (-obj.R0*sind(phi)+ obj.R0*phi*pi/180*cosd(phi))/cosd(obj.beta);
            Y0 =  (obj.R0*cosd(phi)+ obj.R0*phi*pi/180*sind(phi))/cosd(obj.beta);
            [X,Y] = trRot2d(X/cosd(obj.beta),Y+obj.e,X0,Y0,phi);
            drawPolyline(X,Y,'k:')
            drawCross(2*obj.rack.m,2*obj.rack.m,X0,Y0,phi,'m:')
            
        end
        
    end
    
    % calculation methods
    methods (Access = private)
        function [X,Y] = profil(obj,n,t)
            % calculate points on given segment n
            [xi,eta,dydx] = calcProfile(obj.rack,n,t);
            cosb = cosd(obj.beta);
            phi = -(xi + (obj.e + eta)*cosb^2.*dydx)/obj.R0;
            sinp = sin(phi);
            cosp = cos(phi);
            X0 = (obj.R0*phi.*cos(phi) - obj.R0*sin(phi))/cosb;
            Y0 = (obj.R0*cosp + obj.R0*phi.*sinp)/cosb;
            X  = X0 + xi.*cosp/cosb - (eta + obj.e).*sinp;
            Y  = Y0 + xi.*sinp/cosb + (eta + obj.e).*cosp;
        end
        function [X,Y] = profil0(obj,phi)
            % calculate point on involute part
            tana = tand(obj.rack.alpha);
            cosb = cosd(obj.beta);
            sinp = sin(phi);
            cosp = cos(phi);
            xi   = -tana^2*phi*obj.R0/(cosb^2+tana^2) + ...
                (pi*obj.rack.m+4*obj.e*tana)*cosb^2/(4*(cosb^2+tana^2));
            eta  = (pi*obj.rack.m/4 - xi)/tana;
            X0   = (obj.R0*phi.*cosp - obj.R0*sinp)/cosb;
            Y0   = (obj.R0*cosp + obj.R0*phi.*sinp)/cosb;
            X    = X0 + xi.*cosp/cosb - (eta + obj.e).*sinp;
            Y    = Y0 + xi.*sinp/cosb + (eta + obj.e).*cosp;
        end
        function [dX,dY] = profilDer(obj,n,t)
            % calculate derivatives of tooth profile
            narginchk(3,3)
            nargoutchk(2,2)
            %check input
            validateattributes(n, {'numeric'}, {'>',0,'<',5,'integer','scalar'});
            validateattributes(t, {'numeric'}, {'real','vector'});
            % calculate derivatives of rack profile
            [xi,eta,dydx,dx,dy,ddx,ddy] = calcProfile(obj.rack,n,t);
            cosb = cosd(obj.beta);
            phi  = - (xi + (eta + obj.e)*cosb^2.*dydx)/obj.R0;
            dp   = -(((obj.e + eta).*ddy + dy.^2)*cosb^2 + ...
                (obj.R0*phi + xi).*ddx + dx.^2)./(obj.R0*dx);
            cosp = cos(phi);
            sinp = sin(phi);
            dX   = -((obj.e*cosb + eta*cosb).*cosp + (obj.R0*phi + xi).*sinp).*dp/cosb + ...
                cosp.*dx/cosb - sinp.*dy;
            dY   = ((obj.R0*phi + xi).*cosp/cosb - (obj.e + eta).*sinp).*dp + ...
                sinp.*dx/cosb + cosp.*dy;
        end
        function [dX,dY] = calcTangent(obj,n,t)
            % calculate tangent to tooth profile
            narginchk(3,3)
            nargoutchk(2,2)
            %check input
            validateattributes(n, {'numeric'}, {'>',0,'<',5,'integer','scalar'});
            validateattributes(t, {'numeric'}, {'real','vector'});
            % calculate derivatives of rack profile
            [xi,eta,dydx,dx,dy,ddx,ddy] = calcProfile(obj.rack,n,t);
            cosb = cosd(obj.beta);
            phi  = - (xi + (eta + obj.e)*cosb^2.*dydx)/obj.R0;
            dp   = -(((obj.e + eta).*ddy + dy.^2)*cosb^2 + ...
                (obj.R0*phi + xi).*ddx + dx.^2)./(obj.R0*dx);
            cosp = cos(phi);
            sinp = sin(phi);
            dX   = -((obj.e*cosb + eta*cosb).*cosp + (obj.R0*phi + xi).*sinp).*dp/cosb + ...
                cosp.*dx/cosb - sinp.*dy;
            dY   = ((obj.R0*phi + xi).*cosp/cosb - (obj.e + eta).*sinp).*dp + ...
                sinp.*dx/cosb + cosp.*dy;
            % normalize
            aXY = sqrt(dX.^2 + dY.^2);
            dX = -dX./aXY;
            dY = -dY./aXY;
        end
        function val = calcLength(obj)
            % calaclulate length of each profile segment
            val = zeros(4,1);
            for n = 1:4
                if n == 1 && obj.te_(n)  <= 0
                    val(n) = 0;
                    continue
                end
                try
                    %====================================
                    val(n) = integral(@fun,obj.ts_(n),obj.te_(n));
                    %====================================
                catch
                    val(n) = NaN;
                end
            end
            function v = fun(t)
                [dX,dY] = profilDer(obj,n,t);
                v = sqrt(dX.^2 + dY.^2);
            end
        end
        function [ts,te] = calcParam(obj)
            % calculate start and end parameters for tooth profile
            % this can be done analyticaly
            tana = tand(obj.rack.alpha);
            cosb = cosd(obj.beta);
            phib = -tana/cosb - (obj.rack.m*pi/4 + obj.e*tana)/obj.R0;
            phi  = phib + sqrt((obj.Ra/obj.Rb)^2 - 1);
            xi   =(-phi*obj.R0*tana^2 + (obj.rack.m*pi/4 + ...
                obj.e*tana)*cosb^2)/(tana^2 + cosb^2);
            t20  = (obj.rack.m*obj.rack.u*tana - obj.rack.m*pi/4 + xi)/...
                (tana*obj.rack.m*(obj.rack.u+1));
            [X,Y] = profil0(obj,phi);
            t11 = -obj.R0*atan(X/Y)/obj.rack.m/(obj.rack.u*tana - pi/4);
            %=====================================================
            %tt = mfsolve(@fun12,[0.5,0.5]); %,'Display','none');
            %=====================================================
            %t11 = tt(1);
            %t20 = tt(2);
            % If there is intersection between involute and fillet - calculate it
            tc1 = (cosd(obj.beta)*sind(obj.rack.alpha)^2*obj.R0+cosd(obj.beta)*obj.e + ...
                obj.rack.u)/(1 + obj.rack.u);
            if tc1 < 1
                try
                    %=============================
                    zz = mfsolve(@fun23,[t20,0]);
                    %=============================
                    t21 = zz(1);
                    t30 = zz(2);
                catch
                    t21 = 1;
                    t30 = 0;
                end
            else
                t21 = 1;
                t30 = 0;
            end
            t21 = min(1,t21);
            t30 = max(0,t30);
            % start and end parameters of segments
            ts = [ 0,  t20, t30, 0];
            te = [t11, t21,   1, 1];
            %{
            function z = fun12(tt)
                % intersection between seg1 and seg2
                t1 = tt(1);
                t2 = tt(2);
                [tx1,ty1] = profil(obj,1,t1);
                [tx2,ty2] = profil(obj,2,t2);
                z(1) = tx1 - tx2;
                z(2) = ty1 - ty2;
            end
            %}
            function z = fun23(tt)
                % intersection between seg2 and seg3
                t2 = tt(1);
                t3 = tt(2);
                [tx2,ty2] = profil(obj,2,t2);
                [tx3,ty3] = profil(obj,3,t3);
                z(1) = tx2 - tx3;
                z(2) = ty2 - ty3;
            end
        end
        function [X,Y] = calcPoints(obj,nz)
            % cacultae points on gear for nz teeth
            if nargin < 2
                nz = obj.z;
            else
                validateattributes(nz,{'numeric'},{'>',0,'<=',obj.z,'integer','scalar'});
            end
            dsmin = obj.Lt/obj.np;
            [X,Y] = calcPoints1(obj,dsmin,nz);
            %{
            X = [];
            Y = [];
            ts = obj.tpar(:,1);
            te = obj.tpar(:,2);
            for n = 1:4
                nps = max(3,obj.np*obj.L(n)/obj.Lt);
                [XX,YY] = profil(obj,n,linspace(ts(n),te(n),nps));
                X = [X;XX'];
                Y = [Y;YY'];
            end
            X = [-flipud(X);X];
            Y = [ flipud(Y);Y];
            XX = X;
            YY = Y;
            for n = 1:nz - 1
                [X,Y] = trRot2d(X,Y,0,0,360/obj.z);
                XX = [X;XX];
                YY = [Y;YY];
            end
            X = XX;
            Y = YY;
            [X,Y] = deleteDuplicate(X,Y);
            X = flipud(X);
            Y = flipud(Y);
            %}
        end
        function [X,Y] = calcPoints1(obj,dsmin,nz)
            % caculte points on gear for nz teeth  where distance between points is dsmin
            if nargin < 3
                nz = obj.z;
            else
                validateattributes(nz,{'numeric'},{'>',0,'<=',obj.z,'integer','scalar'});
            end
            X = [];
            Y = [];
            % [ts,te] = calcParam(obj);
            for n = 1:4
                % how many points ?
                if obj.L(n) == 0  % Avoid pointed teeth
                    continue
                end
                nps = max(1,ceil(obj.L(n)/dsmin)) + 1;
                if n == 1 || n == 4
                    [XX,YY] = profil(obj,n,linspace(obj.ts_(n),obj.te_(n),nps)');
                    %  scatter(XX,YY,30,'r')
                else
                    XX = zeros(nps,1);
                    YY = zeros(nps,1);
                    ds = obj.L(n)/nps;
                    for k = 1:nps+1
                        if k == 1
                            t = obj.ts_(n);
                        else
                            [dX,dY] = profilDer(obj,n,t);
                            t = t + ds/sqrt(dX^2+ dY^2);
                        end
                        [XX(k),YY(k)] = profil(obj,n,t);
                    end
                end
                X = [X;XX];
                Y = [Y;YY];
            end
            X = [-flipud(X);X];
            Y = [ flipud(Y);Y];
            XX = X;
            YY = Y;
            for n = 1:nz - 1
                [X,Y] = trRot2d(X,Y,0,0,360/obj.z);
                XX = [X;XX];
                YY = [Y;YY];
            end
            X = XX;
            Y = YY;
            [X,Y] = deleteDuplicate(X,Y);
            X = flipud(X);
            Y = flipud(Y);
            % rotate in CW direction
            if nz < obj.z && nz > 2
                [X,Y] = trRot2d(X,Y,0,0,-360/obj.z);
            end             
        end
        function R = calcRadius(obj)
            % calculate radius corespond to the key points
            %[ts,te] = calcParam(obj);
            for n = 1:4
                [X,Y] = profil(obj,n,[obj.ts_(n),obj.te_(n)]);
                R(2*n-1:2*n) = sqrt(X.^2 + Y.^2);
            end
        end

    end
    
    % calculation of tooth thickness
    methods (Access = private)
        function s = toothThickness(obj,R)
            % tooth thickness in involute segment
            tana  = tand(obj.rack.alpha);
            cosb  = cosd(obj.beta);
            phib  = -tana/cosb - 1/obj.R0*(obj.rack.m*pi/4 + obj.e*tana);
            phi   = phib+sqrt((R/obj.Rb)^2 - 1);
            [X,Y] = profil0(obj,phi);
            theta = atan(X/Y);
            s = 2*R*theta;
        end
        function s = toothThickness1(obj,R)
            % numerical calculation of tooth thickness
            validateattributes(R,{'numeric'},{'>=',obj.Rd,'<=',obj.Ra,'real','scalar'});
            s = NaN;
            if R == obj.Ra
                n = 2;
                t = obj.ts_(n);
            elseif R == obj.Rd
                n = 4;
                t = obj.ts_(n);
            elseif R == obj.Ru
                n = 2;
                t = obj.te_(2);
            else
                if R > obj.Ru
                    n = 2;
                    try
                        %==============================
                        t = mfzero(@fun,[obj.ts_(n),obj.te_(n)]);
                        %==============================
                    catch
                        return
                    end
                else
                    n = 3;
                    try
                        %==============================
                        t = mfzero(@fun,[obj.ts_(n),obj.te_(n)]);
                        %===============================
                    catch
                        return
                    end
                end
            end
            [X,Y] = profil(obj,n,t);
            theta = atan(X/Y);
            s = 2*R*theta;
            function val = fun(t)
                [X,Y] =  profil(obj,n,t);
                val = X^2 + Y^2 - R^2;
            end
        end
    end
    
    % var. methods
    methods (Access = private)
        function [x,y] = point2(obj,t)
            % calculate a path of rack point #2
            xx = obj.rack.m*(pi/4 - obj.rack.u*tand(obj.rack.alpha));
            yy = obj.rack.u*obj.rack.m;
            dydx = -t/tand(obj.rack.alpha);
            phi = -(xx + (obj.e + yy)*cosd(obj.beta)^2.*dydx)/obj.R0;
            x   = (-(obj.R0 + (obj.e + yy)*cosd(obj.beta)).*sin(phi) + ...
                (obj.R0*phi + xx).*cos(phi))/cosd(obj.beta);
            y   = ( (obj.R0 + (obj.e + yy)*cosd(obj.beta)).*cos(phi) + ...
                (obj.R0*phi + xx).*sin(phi))/cosd(obj.beta);
        end
    end
    
    methods 
        function [r,tc] = calcRc(obj)
            % calculate critical radius for pointed teeth
            try
                tc = mfsolve(@fun,0);
                [X,r] = profil(obj,2,tc);
            catch
                r  = NaN;
                tc = NaN;
            end
            function val = fun(t)
                [X,~] = profil(obj,2,t);
                val = X;
            end
        end
        function val = calcXmax(obj)
            % max. profile shif coefficient
            tana = tand(obj.rack.alpha);
            cosb = cosd(obj.beta);
            %find interval
            phi0=acot(tana/cosb);
            % solve equation
            try
                %===================================
                phi = mfzero(@fun,[0,0.99*phi0]);
                %==================================
                val = -1/tana*(obj.R0*(phi - sin(phi)*...
                    (cosb^2 + tana^2)/cosb/(cosb*cos(phi) - tana*sin(phi)))...
                    + pi*obj.rack.m/4);
            catch
                val = NaN;
            end
            function f = fun(phi)
                cosp = cos(phi);
                sinp = sin(phi);
                a1 = ((1 - cosp)*tana - cosb*sinp)./(cosb*cosp - tana*sinp);
                a2 = obj.rack.m/obj.R0*(pi/4 - obj.rack.u*tana);
                f = phi + a1 + a2;
            end
        end
    end
    
end

function x = mfsolve(fun,x0)
% Wrapper to fsolve
    tol = 1.e-6;
    x = fsolve(fun,x0); %lsqnonlin
    f = fun(x);
    if norm(f) > tol
        error('Fail to solve equations.')
    end
end

function x = mfzero( fun, x0)
% Wrapper to fzero
    x = fzero(fun,x0);
end

function [dd,mm,ss] = deg2dms(deg)
%DMS2DEG  Convert decimal degrees to ddd:mm:ss
    dd = fix(deg);
    mm = fix((deg - dd)*60); 
    ss = round((deg - dd - mm/60)*3600,2);
    if ss > 59.9999
        mm = mm + 1;
        ss = 0;
    end
    if mm == 60
        dd = dd + 1;
        mm = mm - 60;
    end
end

function [x,y] = deleteDuplicate(X,Y)
%DELETEDUPLICATE Delete duplicate entries in X,Y arrays which forms a
%contour
    tol = 1e-4;
    narginchk(2,2)
    nargoutchk(2,2)
    validateattributes(X, {'numeric'}, {'real','vector'});
    validateattributes(Y, {'numeric'}, {'real','vector'}); 
    if ~isequal(size(X),size(Y))
        error('Arrays must be of the same size.')
    end
    np = length(X);
    x  = nan(size(X));
    y  = nan(size(Y));
    k  = 1;
    x(1) = X(1);
    y(1) = Y(1);
    for n = 2:np 
        if (X(n) - x(k))^2 + (Y(n) - y(k))^2 < tol^2
            continue
        end
        k = k + 1;
        x(k) = X(n);
        y(k) = Y(n);        
    end
    if (x(1) - x(k))^2 + (y(1) - y(k))^2 < tol^2
        k = k - 1;
    end
    x = x(1:k);
    y = y(1:k);
end



