classdef qtn  
    % qtn  Basic constructor for quaternions. 
    %
    %   qtn(a,b,c,d)
    %   qtn(a,[v1 v2 v3])
    %   qtn(a) 
    %
    %   Methods
    %
    %   plus        p + q
    %   minus       p - q
    %   uminus      -q
    %   mtimes      p*q
    %   mldivide    p\q
    %   mrdivide    p/q 
    %   mpower      q^n
    %   inv(q)
    %   norm(q)
    %   disp(q)
    %
    %   Functions
    %
    %   rot
    %   Euler
    %   randq
    %   qtn_examples
    %
    %   Notes:
    %      All angles are in degrees.
    %      This object is not vectorized.
    
    %   Copywrite 2021  Ed Angel and Cleve Moler

    properties
        s  % scalar part, scalar
        v  % vector part, length = 3
    end

    methods
        function q = qtn(a,b,c,d)
            if nargin == 0
                q.s = [];
                q.v = [];
            elseif isa(a,'qtn')
                q = a;
            elseif isa(a,'struct')
                q.s = a.s;
                q.v = a.v;
            elseif nargin == 4
                q.s = a;
                q.v = [b c d];
            elseif nargin == 2
                q.s = a;
                q.v = b;
            elseif isscalar(a)
                q.s = a;
                q.v = [0 0 0];
            else
                q.s = a(1);
                q.v = a(2:end);
            end
        end

        function r = plus(p,q)
            % p+q
            p = qtn(p);
            q = qtn(q);
            r.s = p.s + q.s;
            r.v = p.v + q.v;
            r  = qtn(r);
        end

        function r = minus(p,q)
            % p-q
            p = qtn(p);
            q = qtn(q);
            r.s = p.s - q.s;
            r.v = p.v - q.v;
            r  = qtn(r);
        end
        
        function r = uminus(p)
            % -p
            r.s = -p.s;
            r.v = -p.v;
            r  = qtn(r);
        end
       
        function r = mtimes(p,q)
            % p*q
            p = qtn(p);
            q = qtn(q);
            r.s = p.s*q.s - dot(p.v,q.v);
            r.v = p.s*q.v + q.s*p.v + cross(p.v,q.v);
            r = qtn(r);
        end

        function x = mldivide(p,q)
            % p\q without inv(p)
            p = qtn(p);
            q = qtn(q);
            p.v = -p.v;
            x = mtimes(p,q);
            x = [x.s x.v]/sum([p.s p.v].^2);
            x = qtn(x);
        end

        function x = mrdivide(p,q)
            % p/q without inv(q);
            p = qtn(p);
            q = qtn(q);
            q.v = -q.v;
            x = mtimes(p,q);
            x = [x.s x.v]/sum([q.s q.v].^2);
            x = qtn(x);
        end

        function x = mpower(q,n)
            % q^n
            x = qtn(1);
            if n > 0
                for k = 1:n
                    x = mtimes(q,x);
                end
            else
                for k = 1:(-n)
                    x = mldivide(q,x);
                end
            end
        end
        
        function r = inv(p)
            % inv(p)
            r = [p.s, -p.v];
            r = r/sum(r.^2);
            r = qtn(r);
        end

        function r = norm(p)
            % norm(p)
            r = sqrt(sum([p.s, p.v].^2));
        end

        function disp(q)
            % disp(q)
            q = round([q.s q.v],15);
            pm = blanks(4);
            pm(q >= 0) = '+';
            pm(q < 0) = '-';
            a = abs(q);
            fprintf('   %.6g %1s %.6gi %1s %.6gj %1s %.6gk\n', ...
                q(1),pm(2),a(2),pm(3),a(3),pm(4),a(4))
        end
    end
end