classdef gearsInMesh < handle
    %GEARSINMESH Create gearsInMesh object
    %   
    
    properties
        G1      % gear 1
        G2      % gear 2
    end

    properties (Dependent)
        a      % center distance
        alphaw % working pressure angle in degrees
        cw     % working tip tooth clearance
        Rw1    % working pitch circle radius for gear 1
        Rw2    % working pitch circle radius for gear 2
        invalid
    end
    properties (Access = private)
        a_
        alphaw_
        cw_
        invalid_ % set to true if fail to calculate a_
        Rw1_
        Rw2_
    end
    
        
    methods
        function obj = gearsInMesh(G1,G2)
            if ~(isa(G1,'gear') && isa(G2,'gear'))
                error('Input must be gear objects.')
            end
            if G1.rack.m ~= G2.rack.m || G1.rack.alpha ~= G2.rack.alpha
                error('Gears in mesh must have the same rack.')
            end
            obj.G1 = G1;
            obj.G2 = G2;
            ainit(obj)
        end
        
    end
    
    methods (Access = private)
        function ainit(obj)
            obj.alphaw_ = workingPressureAngle(obj.G1.rack.alpha,...
                obj.G1.z,obj.G2.z,obj.G1.x,obj.G2.x);
            obj.a_ = centerDistance(obj.G1.rack.alpha,obj.G1.rack.m,...
                obj.G1.z,obj.G2.z,obj.G1.x,obj.G2.x,obj.G1.beta);
            obj.cw_ = obj.a_ - obj.G1.Rr - obj.G2.Ra + obj.G1.rack.m*(1 - obj.G1.x) ...
                + obj.G1.rack.c;
            obj.Rw1_ = obj.G1.Rb/cosd(obj.alphaw);
            obj.Rw2_ = obj.G2.Rb/cosd(obj.alphaw);
            if obj.a_ < 0 || obj.a_> obj.G1.Ra + obj.G2.Ra
                warning('Invalid center distance. Set to default\n Correct profile shifts.')
                obj.a_ = obj.G1.Rr + obj.G2.Rr;
                obj.alphaw_ = obj.G1.rack.alpha;
                obj.cw_ = obj.a_ - obj.G1.Rr - obj.G2.Ra + obj.G1.rack.m ...
                + obj.G1.rack.c;
                obj.Rw1_ = obj.G1.Rb/cosd(obj.alphaw);
                obj.Rw2_ = obj.G2.Rb/cosd(obj.alphaw);
                obj.invalid_ = true;
            else
                obj.invalid_ = false;
            end
        end
    end
    
    methods
        function out = get.a(obj)
            out = obj.a_;
        end
        function out = get.alphaw(obj)
            out = obj.alphaw_;
        end
        function out = get.cw(obj)
            out = obj.cw_;
        end
        function out = get.Rw1(obj)
            out = obj.Rw1_;
        end
        function out = get.Rw2(obj)
            out = obj.Rw2_;
        end
        function out = get.invalid(obj)
            out = obj.invalid_;
        end        
    end
    
    methods
        function plot(obj,varargin)
            %Input--
            %   obj  -- gearsInMesh object
            %
            %Optional arguments:
            %   'save'     -- save figure as jpg
            %
            %Optional name-value pairs:
            %   '-a',a     -- center distance                
            %   '-zoom',fc -- scale factor in units of module. 2 means window of size
            %                   2 modules at contact point . Default is no zoom
            %   '-th1',th1 -- rotation angle in degrees for gear 1
            %   '-np',np     -- number of points for gear contour approx.
            
            narginchk(1,inf)
            nargoutchk(0,0)
                       
            ainit(obj);
            
            m = obj.G1.rack.m;
            sz = m/4; % size of point marker
            
            % default values
            fc  = 0;
          %  a  = centerDistance(obj.G1.rack.alpha,obj.G1.rack.m,...
          %      obj.G1.z,obj.G2.z,obj.G1.x,obj.G2.x,obj.G1.beta);
            aa = obj.a_;
            th1 = 0;
            np   = 20;  % # of points to draw the profile
            save = false;
            
            if ~isempty(varargin)
                for k = 1:length(varargin)
                    if ~ischar(varargin{k})
                        continue
                    end
                    switch lower(varargin{k})                    
                        case {'-fc','-zoom'}
                            fc = varargin{k + 1};
                            validateattributes(fc,  {'numeric'}, {'>=',2,'real','scalar'});
                        case '-th1'
                            th1 = varargin{k + 1};
                            validateattributes(th1,{'numeric'}, {'real','scalar'});
                        case '-np'
                            np = varargin{k + 1};
                            validateattributes(np,{'numeric'}, {'positive','integer','scalar'});
                        case '-a'
                            aa = varargin{k + 1};
                            validateattributes(aa,{'numeric'}, {'>',0,'real','scalar'});
                        case 'save'
                            save = true;
                        otherwise
                    end
                end
            end
            th2 = -th1*obj.G1.Rr/obj.G2.Rr;
            
            [X1,Y1] = gearContour(obj.G1,'-np',np);
            [X2,Y2] = gearContour(obj.G2,'-np',np);
            
            [X1,Y1] = trRot2d(X1,Y1,0,0, -90);
            [X2,Y2] = trRot2d(X2,Y2,0,0,  90);
            [X2,Y2] = trRot2d(X2,Y2,0,0,  180/obj.G2.z);
            
            [X1,Y1] = trRot2d(X1,Y1,0,0,  th1);
            [X2,Y2] = trRot2d(X2,Y2,0,0,  th2);
            
            drawInit;
            set(gcf, 'Position', [100 100 600 600]);
            drawSet('LineWidth',2,'LineColor','k');
            drawPolyline(X2+aa,Y2)
            drawPolyline(X1,Y1)
            drawPoint(1,sz,0,0)
            drawPoint(1,sz,aa,0)
          %  title(sprintf('m = %g, z_1 = %g, z_2 = %g, a = %g',m,obj.G1.z,obj.G2.z,aa),...
          %      'FontSize',12,'FontWeight','normal')  
                title(sprintf('m = %g, z_1 = %g, z_2 = %g, x_1 + x_2 = %g, a = %g',...
                    m,obj.G1.z,obj.G2.z,obj.G1.x+obj.G2.x,aa),...
                    'FontSize',12,'FontWeight','normal')            
            if fc > 0
                xmin = -fc*m+obj.G1.Rr;
                xmax = fc*m + obj.G1.Rr;
                ymin = -fc*m;
                ymax = fc*m;
                drawLimits(xmin,xmax,ymin,ymax)
            end
            axis off
            if save
                fn = drawSave;
                fprintf('Figure is saved in file %s\n.',fn)
            end            
        end
    end

    methods
        function animate(obj,varargin)
            %ANIMATE Draw gears in mesh in successive positions
            %
            %Input--
            %   obj  -- gearsInMesh object
            %
            %Optional name-value pairs:
            %   '-a',a     -- center distance
            %   '-zoom',fc -- scale factor in units of module. 2 means window of size
            %                   2 modules at contact point . Default is no zoom
            %   '-dth1',dth1 -- gear 1 rotation step size in degrees
            %   '-nr',nr     -- number of turns before stop
            %   '-step','on' -- interactive step using arrow or mouse. THIS DOES NOT
            %                   WORK WELL !!!
            
            narginchk(1,inf)
            nargoutchk(0,0)
               
            ainit(obj);
            
            m = obj.G1.rack.m;
            sz = m/4; % size of point marker
            
            % default values
            fc   = 0;
            %a  = centerDistance(obj.G1.rack.alpha,obj.G1.rack.m,...
            %    obj.G1.z,obj.G2.z,obj.G1.x,obj.G2.x,obj.G1.beta);
            aa = obj.a_;
            dth1 = 1;
            nr   = 1;  % # of turns
            step = false;
            save = false;
            
            if ~isempty(varargin)
                for k = 1:length(varargin)
                    switch lower(varargin{k})
                        case {'-fc','-zoom'}
                            fc = varargin{k + 1};
                            validateattributes(fc,  {'numeric'}, {'>=',2,'real','scalar'});
                        case '-dth1'
                            dth1 = varargin{k + 1};
                            validateattributes(dth1,{'numeric'}, {'real','scalar'});
                        case '-nr'
                            nr = varargin{k + 1};
                            validateattributes(nr,{'numeric'}, {'positive','real','scalar'});
                        case '-a'
                            aa = varargin{k + 1};
                            validateattributes(aa,{'numeric'}, {'>',0,'real','scalar'});
                        case '-step'
                            % this does not work well !!!!!
                            if isequal(lower(varargin{k + 1}),'on')
                                step = true;
                                save = false;
                            end
                        case 'save'
                            save = true;
                            step = false;
                        otherwise
                    end
                end
            end
            nn = abs(nr*360/dth1);
            dth2 = -dth1*obj.G1.Rr/obj.G2.Rr;
            
            
            [X1,Y1] = gearContour(obj.G1,'-nn',20);
            [X2,Y2] = gearContour(obj.G2,'-nn',20);
            
            [X1,Y1] = trRot2d(X1,Y1,0,0, -90);
            [X2,Y2] = trRot2d(X2,Y2,0,0,  90);
            [X2,Y2] = trRot2d(X2,Y2,0,0,  180/obj.G2.z);
            
            fg = drawInit;
            set(gcf, 'Position', [100 100 600 600]);
                        drawSet('LineWidth',2,'LineColor','k');
            
            if save
                fprintf('Start animatin file...\n')
                v = VideoWriter('gearsInMesh.avi');
                open(v)
            end
            for n = 1:nn
                [X1,Y1] = trRot2d(X1,Y1,0,0,  dth1);
                [X2,Y2] = trRot2d(X2,Y2,0,0,  dth2);
                drawInit(fg)
                drawPolyline(X2+aa,Y2)
                drawPolyline(X1,Y1)
                drawPoint(1,sz,0,0)
                drawPoint(1,sz,aa,0)
                if fc > 0
                    xmin = -fc*m+obj.G1.Rr;
                    xmax = fc*m + obj.G1.Rr;
                    ymin = -fc*m;
                    ymax = fc*m;
                    drawLimits(xmin,xmax,ymin,ymax)
                end
                axis off
                title(sprintf('m = %g, z_1 = %g, z_2 = %g, x_1 + x_2 = %g, a = %g\n %d/%d',...
                    m,obj.G1.z,obj.G2.z,obj.G1.x+obj.G2.x,aa,n,nn),...
                    'FontSize',12,'FontWeight','normal')
                if save
                       frame = getframe(gcf);
                        writeVideo(v,frame);
                        continue
                end
                if step
                    title('use arrows, esc to quit')
                    [~,~,key] = ginput(1);
                    switch key
                        case {1,29,30} % arrow up/ right
                            if dth1 < 0
                                dth1 = -dth1;
                                dth2 = -dth2;
                            end
                        case {2,27} % esc
                            step = false;
                        case {3,28,31} % arrow down/left
                            if dth1 > 0
                                dth1 = -dth1;
                                dth2 = -dth2;
                            end
                        otherwise
                    end
                end
            end
            if save
                close(v);
                delete(v); % !!!!
                fprintf('Video svaed in %s\n','gearsInMesh.avi');
            end
            
        end
    end
    
    methods
        function  print(obj,fid)
            if nargin < 3
                fid = 1;
            end
            ainit(obj);
            fprintf(fid,'\nGears in mesh\n');
            %a  = centerDistance(obj.G1.rack.alpha,obj.G1.rack.m,...
            %    obj.G1.z,obj.G2.z,obj.G1.x,obj.G2.x,obj.G1.beta);
            %aw = workingPressureAngle(obj.G1.rack.alpha,...
            %    obj.G1.z,obj.G2.z,obj.G1.x,obj.G2.x);
            fprintf(fid,'                               Number of teeth:%12d%12d\n',obj.G1.z,obj.G2.z);
            fprintf(fid,'                                        Module:%18.4f       mm\n',obj.G1.rack.m);
            fprintf(fid,'                       Standard pressure angle:%18.4f       \n',obj.G1.rack.alpha);
            fprintf(fid,'                           Tip tooth clearance:%18.4f       mm\n',obj.G1.rack.c);          
            fprintf(fid,'                               Center distance:%18.4f       \n',obj.a);  
            fprintf(fid,'                    Tip shortening coefficient:%12.4f%12.4f\n',obj.G1.u,obj.G2.u);            
            fprintf(fid,'                         Workin pressure angle:%18.4f       \n',obj.alphaw);
            fprintf(fid,'         Sum of the profile shift coefficients:%18.4f       \n',obj.G1.x+obj.G2.x);
            fprintf(fid,'                    Profile shift coefficients:%12.4f%12.4f\n',obj.G1.x,obj.G2.x);            
            fprintf(fid,'Calculated parameters\n');
            fprintf(fid,'               Reference pitch circle diameter:%12.4f%12.4f mm\n',2*obj.G1.Rr,2*obj.G2.Rr);
            fprintf(fid,'                 Working pitch circle diameter:%12.4f%12.4f mm\n',2*obj.Rw1,2*obj.Rw2);
            fprintf(fid,'                      Addendum circle diameter:%12.4f%12.4f mm\n',2*obj.G1.Ra,2*obj.G2.Ra);
            fprintf(fid,'                      Dedendum circle diameter:%12.4f%12.4f mm\n',2*obj.G1.Rd,2*obj.G2.Rd);
            fprintf(fid,'                          Base circle diameter:%12.4f%12.4f mm\n',2*obj.G1.Rb,2*obj.G2.Rb);
            fprintf(fid,'                   Working tip tooth clearance:%18.4f       \n',obj.cw_);            
            fprintf(fid,' Tooth thickness at the reference pitch circle:%12.4f%12.4f mm\n',obj.G1.sr,obj.G2.sr);
            fprintf(fid,'            Tooth thickness at the base circle:%12.4f%12.4f mm\n',obj.G1.sb,obj.G2.sb);
            fprintf(fid,'                           Tip tooth thickness:%12.4f%12.4f mm\n',obj.G1.sa,obj.G2.sa);
            if obj.invalid_
                fprintf('***Warning: ceneter distance and working pressure angle are invalid\n');
                fprintf('   Adjust profile shifts\n')
            end
            if obj.G1.Ru > obj.G1.Rb && obj.G1.z < obj.G1.zmin
                fprintf(fid,'***Warning: gear 1 is undercut.\n');
                fprintf(fid,'   Set profil shift >= %g\n',obj.G1.xmin);
            end
            if obj.G2.Ru > obj.G2.Rb && obj.G2.z < obj.G2.zmin
                fprintf(fid,'***Warning: gear 2 is undercut.\n');
                fprintf(fid,'   Set profil shift >= %g\n',obj.G2.xmin);
            end
            if obj.G1.sa < 0.25*obj.G1.sr
                fprintf(fid,'***Warning: gear 1 has pointed teeth.\n');
            end            
            if obj.G2.sa < 0.25*obj.G2.sr
                fprintf(fid,'***Warning: gear 2 has pointed teeth.\n');
            end
        end
    end
    
end

function a = centerDistance(alpha,m,z1,z2,x1,x2,beta)
%invaw = 2*tand(alpha)*((x1 + x2)/(z1 + z2)) + invd(alpha);
if x1 == x2 || x1 == 0
    a = m*(z1 + z2)/cosd(beta)/2;
    return
end
alphaw = workingPressureAngle( alpha,z1,z2,x1,x2);
y = (z1 + z2)/2*(cosd(alpha)/cosd(alphaw) - 1);
a = ((z1 + z2)/2 + y)*m/cosd(beta);
end

function alphaw = workingPressureAngle(alpha,z1,z2,x1,x2)
if x1 == x2 || x1 == 0 || x1 + x2 == 0
    alphaw = alpha;
    return
end
invaw = 2*tand(alpha)*((x1 + x2)/(z1 + z2)) + invd(alpha);
alphaw = ainvd(invaw);
end

function val = invd(alpha)
    val = tand(alpha) - alpha*pi/180;
end

function val = ainvd(inva)
    a = 0:0.5:89;
    f = fun(a);
    f1 = f(1);
    for n = 2:length(f)
        f2 = f(n);
        if f2*f1 < 0
            val0 = n;
            break
        end
    end
    try    
        val = fzero(@fun,val0);
    catch
        val = NaN;
    end
    function y = fun(x)
        y = invd(x) - inva;
    end
end


