added test_particle

+anakin/basis.m

@@ -34,7 +34,7 @@
   with respect to another basis (symbolic variables must be used)
 * subs: takes values of the symbolic unknowns and returns a basis with
   purely numeric matrix (symbolic variables must be used)    
-* isunitary, isrighthanded: checks the corresponding property and
+* isorthogonal, isrighthanded: checks the corresponding property and
   returns true or false  
 * plot: plots the basis with quiver, at a chosen position
 
@@ -272,13 +272,13 @@ purely numeric matrix (symbolic variables must be used)
         end
     end
     methods % logical tests
-        function isunitary = isunitary(B) % all vectors are unitary and mutually orthogonal
+        function isorthogonal = isorthogonal(B) % all vectors are unitary and mutually orthogonal
             if isa(B.m,'sym') % symbolic inputs
-                isunitary = isAlways(B.m' * B.m == eye(3),'Unknown','false'); % In case of doubt, false
+                isorthogonal = isAlways(B.m' * B.m == eye(3),'Unknown','false'); % In case of doubt, false
             else % numeric input            
-                isunitary = (abs(B.m' * B.m - eye(3))<eps(max(abs(B.m(:))))); 
+                isorthogonal = (abs(B.m' * B.m - eye(3))<eps(max(abs(B.m(:))))); 
             end 
-            isunitary = all(isunitary(:));
+            isorthogonal = all(isorthogonal(:));
         end    
         function isrighthanded = isrighthanded(B) % basis is righthanded
             isrighthanded = (det(B.m) > 0);

+anakin/frame.m

@@ -145,7 +145,7 @@ purely numeric origin point and basis matrix (symbolic variables must be used)
             S_.c = double(subs(S.c,variables,values));
             S_.m = double(subs(S.m,variables,values));
         end         
-    end
+    end 
     methods % plotting
         function h = plot(S,varargin) % plot  
             cc = S.c;
@@ -161,13 +161,16 @@ purely numeric origin point and basis matrix (symbolic variables must be used)
     end
     methods % removed methods 
         function mtimes(~,~)
-            error('This usage of point is not permitted');
+            error('This usage of frame is not permitted');
         end 
         function mrdivide(~,~)
-            error('This usage of point is not permitted');
+            error('This usage of frame is not permitted');
         end 
         function mldivide(~,~)
-            error('This usage of point is not permitted');
+            error('This usage of frame is not permitted');
         end 
+        function isunitary(~) % all vectors are unitary and mutually orthogonal
+            error('This usage of frame is not permitted');
+        end    
     end
 end

+anakin/particle.m

@@ -2,30 +2,32 @@
 particle: class to model a point particle. Inherits from point.
 
 P0 = anakin.particle();  % no arguments return default object  
-P  = anakin.particle(<P|A|a|c|x,y,z>,<mass>,<S1>);
+P  = anakin.particle(<mass>,<P|A|a|c|x,y,z>,<S1>);
 
 where:
 - <> denotes optional arguments
 - | denotes alternative arguments
 - P0 is the default particle (mass 1, located at the origin)
 - P  is a particle 
+- mass is the object mass
 - A is a point
 - a is a vector
 - c is an array with the three vector components
 - x,y,z are the three vector components
-- mass is the object mass
 - S1 is a frame. If given, all previous input as relative to that frame
 
 METHODS: 
 * p: linear momentum in a given reference frame
 * H: angular momentum about a point in a given reference frame
 * T: kinetic energy in a given reference frame
 * inertia: tensor of inertia about a point
+* subs: takes values of the symbolic unknowns and returns a particle
+  object which is purely numeric
 
 AUTHOR: Mario Merino <[email protected]>
 %}
 classdef particle < anakin.point 
-    properties
+    properties (SetAccess = protected)
         mass = 1; % mass of the object
     end
     methods % creation
@@ -39,7 +41,7 @@
                 case 0 % no arguments
                     return;
                 case 1  
-                    Pt = anakin.particle(varargin{1},anakin.frame).c;
+                    Pt = anakin.particle(varargin{1},anakin.frame);
                     P.c = Pt.c; 
                     P.mass = Pt.mass;
                 case 2  
@@ -51,34 +53,32 @@
                             P.mass = varargin{1};
                         end
                     else
-                        Pt = anakin.particle(varargin{1},varargin{2},anakin.frame).c;
+                        Pt = anakin.particle(varargin{1},varargin{2},anakin.frame);
                         P.c = Pt.c;
                         P.mass = Pt.mass;
                     end 
                 case 3  
                     if isa(varargin{end},'anakin.frame') % last is frame
-                        if isa(varargin{1},'anakin.vector') || numel(varargin{1}) == 3 % (vector or components), mass, frame
-                            P.c = varargin{3}.matrix * anakin.vector(varargin{1}).c + varargin{3}.c;
-                            P.mass = varargin{2}; 
+                        if isa(varargin{2},'anakin.vector') || numel(varargin{2}) == 3 % (vector or components), mass, frame
+                            P.c = varargin{3}.matrix * anakin.vector(varargin{2}).c + varargin{3}.c;
+                            P.mass = varargin{1}; 
                         end
                     else
-                        Pt = anakin.particle(varargin{1},varargin{2},varargin{3},anakin.frame).c;
+                        Pt = anakin.particle(varargin{1},varargin{2},varargin{3},anakin.frame);
                         P.c = Pt.c;
                         P.mass = Pt.mass;
                     end  
                 case 4 
-                    if isa(varargin{end},'anakin.frame') % last is frame
-                        if isa(varargin{1},'anakin.vector') || numel(varargin{1}) == 3 % xyz, frame
-                            P.c = varargin{4}.matrix * anakin.vector(varargin{1},varargin{2},varargin{3}).c + varargin{4}.c; 
-                        end
+                    if isa(varargin{end},'anakin.frame') % last is frame 
+                        P.c = varargin{4}.matrix * anakin.vector(varargin{1},varargin{2},varargin{3}).c + varargin{4}.c; 
                     else
-                        Pt = anakin.particle(varargin{1},varargin{2},varargin{3},varargin{4},anakin.frame).c;
+                        Pt = anakin.particle(varargin{1},varargin{2},varargin{3},varargin{4},anakin.frame);
                         P.c = Pt.c;
                         P.mass = Pt.mass;
                     end   
                 case 5 % xyz, mass, frame
-                    P.c = varargin{5}.matrix * anakin.vector(varargin{1},varargin{2},varargin{3}).c + varargin{5}.c;
-                    P.mass = varargin{4}; 
+                    P.c = varargin{5}.matrix * anakin.vector(varargin{2},varargin{3},varargin{4}).c + varargin{5}.c;
+                    P.mass = varargin{1}; 
                 otherwise % other possibilities are not allowed
                     error('Wrong number of arguments in particle');
             end       
@@ -111,6 +111,14 @@
                 S1 = anakin.frame; % default frame
             end
             T = (P.mass/2) * dot(P.vel(S1),P.vel(S1));
+            if isa(T,'sym') % symbolic input
+                T = formula(simplify(T)); % simplify and force sym rather than symfun to allow indexing
+            end
+        end
+        function P_ = subs(P,variables,values) % particularize symbolic vector
+            P_ = P;
+            P_.c = double(subs(P.c,variables,values));
+            P_.mass = double(subs(P.mass,variables,values));
         end
     end      
 end

tests/test_particle.m

@@ -0,0 +1,158 @@
+%{
+Compendium of all tests for the function/class in the name of this file.
+
+NOTES:
+* The functions to be tested must be in the Matlab path. You can run the
+  tests by executing 'runtests'. For convenience, you can set the run
+  botton of Matlab's interface to runtests(<this function>)
+* Your working directory must be the
+  directory where this test file is contained.  
+
+AUTHOR: Mario Merino <[email protected]>
+%}
+function tests = test_particle
+    clc;
+    tests = functiontests(localfunctions);     
+end
+
+%----------------------------------------------------------------------
+%----------------------------------------------------------------------
+%----------------------------------------------------------------------
+
+function test_creator(~) % Call creator without arguments
+    import anakin.*
+    S1 = frame([4,5,6],[0,1,0;-1,0,0;0,0,1]); 
+    x = 1; y = 2; z = 3; c = [x;y;z]; cp = [6;4;9];
+    a = vector(c); A = point(c); 
+    mass = 3;
+
+    P = particle();
+    assert(all(P.pos.components == [0;0;0]));
+    assert(all(P.mass == 1));
+
+    P = particle(particle);
+    assert(all(P.pos.components == [0;0;0]));
+    assert(all(P.mass == 1));
+
+    P = particle(A);
+    assert(all(P.pos.components == c));
+    assert(all(P.mass == 1));
+
+    P = particle(a);
+    assert(all(P.pos.components == c));
+    assert(all(P.mass == 1));
+
+    P = particle(c);
+    assert(all(P.pos.components == c));
+    assert(all(P.mass == 1));
+
+    P = particle(c');
+    assert(all(P.pos.components == c));
+    assert(all(P.mass == 1));
+
+    P = particle(x,y,z);
+    assert(all(P.pos.components == c));
+    assert(all(P.mass == 1));
+
+    P = particle(particle,S1);
+    assert(all(P.pos.components == [4;5;6]));
+    assert(all(P.mass == 1));
+
+    P = particle(A,S1);
+    assert(all(P.pos.components == cp));
+    assert(all(P.mass == 1));
+
+    P = particle(a,S1);
+    assert(all(P.pos.components == cp));
+    assert(all(P.mass == 1));
+
+    P = particle(c,S1);
+    assert(all(P.pos.components == cp));
+    assert(all(P.mass == 1));
+
+    P = particle(c',S1);
+    assert(all(P.pos.components == cp));
+    assert(all(P.mass == 1));
+
+    P = particle(x,y,z,S1);
+    assert(all(P.pos.components == cp));
+    assert(all(P.mass == 1));
+
+    P = particle(mass);
+    assert(all(P.pos.components == [0;0;0]));
+    assert(all(P.mass == mass));
+
+    P = particle(mass,A);
+    assert(all(P.pos.components == c));
+    assert(all(P.mass == mass));
+
+    P = particle(mass,a);
+    assert(all(P.pos.components == c));
+    assert(all(P.mass == mass));
+
+    P = particle(mass,c);
+    assert(all(P.pos.components == c));
+    assert(all(P.mass == mass));
+
+    P = particle(mass,c');
+    assert(all(P.pos.components == c));
+    assert(all(P.mass == mass));
+
+    P = particle(mass,x,y,z);
+    assert(all(P.pos.components == c));
+    assert(all(P.mass == mass));
+
+    P = particle(mass,A,S1);
+    assert(all(P.pos.components == cp));
+    assert(all(P.mass == mass));
+
+    P = particle(mass,a,S1);
+    assert(all(P.pos.components == cp));
+    assert(all(P.mass == mass));
+
+    P = particle(mass,c,S1);
+    assert(all(P.pos.components == cp));
+    assert(all(P.mass == mass));
+
+    P = particle(mass,c',S1);
+    assert(all(P.pos.components == cp));
+    assert(all(P.mass == mass));
+
+    P = particle(mass,x,y,z,S1);
+    assert(all(P.pos.components == cp));
+    assert(all(P.mass == mass));
+
+end
+
+function test_momentum(~) % call momentum, linear momentum, energy
+    import anakin.*
+    if license('test','symbolic_toolbox')
+        syms t;
+        syms theta(t) xi(t);
+        assume([in(t, 'real'), in(theta(t), 'real'), in(xi(t), 'real'),...
+                               in(diff(theta(t),t), 'real'), in(diff(xi(t),t), 'real')]);    
+        c = formula([cos(theta);xi^2;xi]); 
+        cp = formula([-sin(theta)*diff(theta,t);2*xi*diff(xi,t);diff(xi,t)]); 
+        mass = 3;    
+        P = particle(mass,c);
+
+        assert(all(isAlways(P.p.components == 3*cp)));
+        assert(all(isAlways(P.H.components == 3*cross(c,cp))));
+        assert(all(isAlways(P.T == 3*dot(cp,cp)/2)));
+    end  
+end
+
+function test_subs(~) % Particularize a symbolic vector
+    import anakin.*
+    if license('test','symbolic_toolbox') 
+        syms t;
+        syms theta(t) xi(t);
+        assume([in(t, 'real'), in(theta(t), 'real'), in(xi(t), 'real'),...
+                               in(diff(theta(t),t), 'real'), in(diff(xi(t),t), 'real')]);    
+        c = formula([cos(theta);xi^2;xi]);  
+        mass = 3;    
+        P = particle(mass,c);
+
+        assert(P.subs({t,theta,xi},{1,3,-2}).pos == vector([cos(3);4;-2])); % call with cell arrays 
+    end
+end

tests/test_vector.m

@@ -19,35 +19,40 @@
 %----------------------------------------------------------------------
 %----------------------------------------------------------------------
 
-function test_creator(~) % Call creator without arguments
-    import anakin.*
-    B1 = basis([0,1,0],[-1,0,0],[0,0,1]); 
-
-    a = vector; % null vector
-    a = vector(vector); % itself
-    a = vector([1;2;0]); % components in canonical vector basis, column array form
-    a = vector([1,2,0]); % components in canonical vector basis, row array form    
-    a = vector(1,2,0); % components in canonical vector basis, independently given
-    a = vector(vector,B1); % relative vector in another basis
-    a = vector([1;2;0],B1); % components in another basis, column array form
-    a = vector([1,2,0],B1); % components in another basis, row array form       
-    a = vector(1,2,0,B1); % components in another basis, independently given
+function test_creator(~) % Call creator 
+    import anakin.* 
+    B1 = basis([0,1,0],[-1,0,0],[0,0,1]);
+    c = [1;2;0];
+    cp = [-2;1;0];
+
+    assert(vector == vector([0;0;0])); % null vector
+    assert(vector(vector) == vector); % itself
+    assert(all(vector(c).components == c)); % components in canonical vector basis, column array form
+    assert(all(vector(c').components == c)); % components in canonical vector basis, row array form
+    assert(all(vector(c(1),c(2),c(3)).components == c)); % components in canonical vector basis, independently given
+    assert(vector(vector,B1) == vector); % relative vector in another basis
+    assert(vector(c,B1) == vector(cp)); % components in another basis, column array form
+    assert(vector(c',B1) == vector(cp)); % components in another basis, row array form
+    assert(vector(c(1),c(2),c(3),B1) == vector(cp)); % components in another basis, independently given
+end
 
+function test_creator2(~) % Call creator with sym arguments
     if license('test','symbolic_toolbox') 
+        import anakin.*
         syms t theta(t) phi(t);
         assume([in(t, 'real'), in(theta(t), 'real'), in(phi(t), 'real')]);
         B1 = basis([1,0,0],[0,cos(phi),sin(phi)],[0,-sin(phi),cos(phi)]);
-
-        a = vector([cos(theta);sin(theta);0]); % components in canonical vector basis, column array form
-        a = vector([cos(theta),sin(theta),0]); % components in canonical vector basis, row array form
-        a = vector([cos(theta);sin(theta);0],B1); % components in another basis, column array form
-        a = vector([cos(theta),sin(theta),0],B1); % components in another basis, row array form
-        a = vector(cos(theta),sin(theta),0); % components in canonical vector basis, independently given
-        a = vector(cos(theta),sin(theta),0,B1); % components in another basis, independently given
+        c = formula([cos(theta);sin(theta);0]);
+
+        a = vector(c); % components in canonical vector basis, column array form
+        a = vector(c'); % components in canonical vector basis, row array form
+        a = vector(c,B1); % components in another basis, column array form
+        a = vector(c',B1); % components in another basis, row array form
+        a = vector(c(1),c(2),c(3)); % components in canonical vector basis, independently given
+        a = vector(c(1),c(2),c(3),B1); % components in another basis, independently given
     end
 end
 
-
 function test_overloads(~) % components, x,y,z
     import anakin.*