forces and eom

+anakin/frame.m

@@ -80,7 +80,7 @@ purely numeric origin point and basis matrix (symbolic variables must be used)
         function disp(S) % display
             disp('Frame with origin with canonical position:')
             disp(S.origin.pos.components)
-            disp('and basis with canonical rotation matrix:')
+            disp('And basis with canonical rotation matrix:')
             disp(S.basis.matrix)
         end
     end

+anakin/particle.m

@@ -17,12 +17,15 @@
 - S1 is a frame. If given, all previous input as relative to that frame
 
 PROPERTIES:
-* mass
+* mass: the mass of the particle
+* forces: a cell array with all the vector forces acting on the particle
 
 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
+* equations: returns a vector of (symbolic) equations, m*a = F, projected
+  along the vectors of one basis
 * inertia: tensor of inertia about a point
 * subs: takes values of the symbolic unknowns and returns a particle
   object which is purely numeric
@@ -33,6 +36,7 @@
 classdef particle < anakin.point 
     properties  
         mass anakin.tensor = anakin.tensor(1); % mass of the object
+        forces cell = {}; % cell array with all forces acting on the particle
     end
     methods % creation
         function P = particle(varargin) % constructor 
@@ -63,7 +67,7 @@
                     error('Wrong number of arguments in particle');
             end       
         end 
-        function P = set.mass(P,value) % on setting c
+        function P = set.mass(P,value) % on setting mass
             P.mass = value;
             if isa(P.mass,'sym') % symbolic input
                 P.mass = formula(simplify(P.mass)); % simplify and force sym rather than symfun to allow indexing
@@ -72,12 +76,23 @@
                 P.mass = double(P.mass);
             end
         end
+        function P = set.forces(P,value) % on setting forces
+            if ~iscell(value) % validate input
+                value = {value};
+            end
+            for i=1:length(value)
+                if ~isa(value{i},'anakin.tensor') 
+                    error('The forces must be supplied in a cell array of anakin.tensor vectors');
+                end
+            end
+            P.forces = value; 
+        end
     end 
-    methods
+    methods(Hidden = true) % overloads
         function disp(P) % display
             disp('Particle with mass:')
             disp(P.mass.components)            
-            disp('canonical position vector components:')
+            disp('And canonical position vector components:')
             disp(P.pos.components)            
         end
     end
@@ -107,6 +122,24 @@ function disp(P) % display
             end
             T = (P.mass/2) * dot(vel,vel); 
         end
+        function eqs = equations(P,B1) % returns vector of equations of motion projected in basis B1
+            MA = P.mass*P.accel;
+            F = anakin.tensor([0;0;0]); % allocate;
+            for i=1:length(P.forces)
+                F = F + P.forces{i};
+            end
+            if ~exist('B1','var')
+                B1 = anakin.basis;
+            elseif isa(B1,'anakin.frame') 
+                B1 = B1.basis; % extract basis
+            end
+            eqs = sym(zeros(B1.spacedim,1));
+            for i=1:B1.spacedim
+                ma = MA * B1.e(i);
+                f = F * B1.e(i);
+                eqs(i) = (ma.components == f.components);
+            end                
+        end
         function inertia = inertia(P,O) % inertia tensor of the particle with respect to point O
             if exist('O','var')
                 r = P.pos0 - O;

+anakin/point.m

@@ -110,10 +110,12 @@ function disp(A) % display
     end 
     methods % plotting
         function h = plot(A,varargin) % plot
-            cc = A.pos0.components;
-            hold on
-            h = line(cc(1),cc(2),cc(3),'color','k','marker','.','linestyle','none');
-            hold off
+            if A.spacedim ~= 3
+                error('This functionality is only available for points in 3D space');
+            end
+            c = A.pos0.components; 
+            h = line;
+            set(h,'XData',c(1),'YData',c(2),'ZData',c(3),'color','k','marker','.','linestyle','none');
             if ~isempty(varargin)
                 set(h,varargin{:}); % set options stored in varargin
             end