quaternion.m

classdef quaternion
    %quaternion Summary of this class goes here
    %   Detailed explanation goes here
    
    properties
        coefficients
        bases
        expression
        MultTable
    end
    
    methods
        function obj = quaternion(base)
            %quaternion Construct an instance of this class
            %   Detailed explanation goes here
            assert(isstring(base))
            q0 = sym(base+"_0");
            q1 = sym(base+"_1");
            q2 = sym(base+"_2");
            q3 = sym(base+"_3");
            obj.coefficients = [q0 q1 q2 q3];
            
            
            e0 = sym("e_0");
            e1 = sym("e_1");
            e2 = sym("e_2");
            e3 = sym("e_3");
            obj.bases = [e0 e1 e2 e3];
            
            MT = [e0  e1  e2  e3;...
                e1 -e0  e3 -e2;...
                e2 -e3 -e0  e1;...
                e3  e2 -e1 -e0];
            
            obj.MultTable = MT;
            obj = obj.buildExpression();
        end
        
        function obj = product(obj,Q1)
            %METHOD1 Summary of this method goes here
            %   Detailed explanation goes here
            assert(strcmpi(class(Q1), class(obj)))
            prod = sym(0);
            for n1 = 1:4
                for n2 = 1:4
                    basis = obj.MultTable(n1,n2);
                    idx = (obj.coefficients == basis);
                    prod = prod + obj.coefficients(n1) * Q1.coefficients(n2) * basis;
                end
            end
            obj.expression = prod;
            obj = obj.expression_ExtractCoefficients();
        end
        
        function obj = buildExpression(obj)
            obj.expression = obj.coefficients .* obj.bases;
        end
        
        function obj = expression_ExtractCoefficients(obj)
            % Reset coefficients to 0
            obj.coefficients = sym(zeros(1,4));
            
            collectTerms = collect(obj.expression,obj.bases);
            formParts = children(collectTerms);
            if any(ismember(formParts,obj.bases))
                %%% Quaternion only in R^1 %%%
                % Get the index of formParts that is the basis
                fIDX = ismember(formParts,obj.bases);
                bIDX = ismember(obj.bases,symvar(formParts(fIDX)));
                formParts = formParts(~fIDX);
                formParts = simplify(formParts);
                obj.coefficients(bIDX) = formParts;
            else
                missingBases = true(1,length(obj.bases));
                for ii = 1:4
                    if ii > length(formParts)
                        % Not all bases are present, i.e. multiplied by 0
                    else
                        bIDX = ismember(obj.bases,symvar(formParts(ii)));
                        missingBases(bIDX) = false;
                        
                        formParts(ii) = subs(formParts(ii),obj.bases(bIDX),sym(1));
                        formParts(ii) = simplify(formParts(ii));
                        obj.coefficients(bIDX) = formParts(ii);
                    end
                end
                obj.coefficients(missingBases) = sym(0);
            end
            
            obj = obj.buildExpression();
        end
        
        function obj = conjugate(obj)
            obj.coefficients(2:end) = -obj.coefficients(2:end);
            obj = obj.buildExpression();
        end
    end
end