Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[8770] trunk/octave-forge/main/mechanics/inst
|
From: <jpi...@us...> - 2011-10-18 06:38:15
|
Revision: 8770
http://octave.svn.sourceforge.net/octave/?rev=8770&view=rev
Author: jpicarbajal
Date: 2011-10-18 06:38:04 +0000 (Tue, 18 Oct 2011)
Log Message:
-----------
mechanics. Renewing interfaace. Deprecating some functions
Modified Paths:
--------------
trunk/octave-forge/main/mechanics/inst/@rigidbody/plot.m
trunk/octave-forge/main/mechanics/inst/@rigidbody/rigidbody.m
trunk/octave-forge/main/mechanics/inst/principalaxes.m
Added Paths:
-----------
trunk/octave-forge/main/mechanics/inst/deprecated/
trunk/octave-forge/main/mechanics/inst/deprecated/area_poly2d.m
trunk/octave-forge/main/mechanics/inst/deprecated/center_mass_poly2d.m
trunk/octave-forge/main/mechanics/inst/deprecated/inertia_moment_ncpoly2d.m
trunk/octave-forge/main/mechanics/inst/deprecated/inertia_moment_poly2d.m
trunk/octave-forge/main/mechanics/inst/deprecated/second_moment_poly2d.m
trunk/octave-forge/main/mechanics/inst/inertiamoment.m
trunk/octave-forge/main/mechanics/inst/masscenter.m
Modified: trunk/octave-forge/main/mechanics/inst/@rigidbody/plot.m
===================================================================
--- trunk/octave-forge/main/mechanics/inst/@rigidbody/plot.m 2011-10-17 20:26:22 UTC (rev 8769)
+++ trunk/octave-forge/main/mechanics/inst/@rigidbody/plot.m 2011-10-18 06:38:04 UTC (rev 8770)
@@ -27,8 +27,9 @@
if has_shape
shapedata = obj.Shape.Data;
if iscell (shapedata)
- error('rigidBody.plot:Devel','Polyline shape, not yet implemented.');
+ shapedata = shape2polygon (shapedata);
end
+
# Rotate
R = cos(obj.Angle)*eye(2) + sin(obj.Angle)*[0 1; -1 0];
shapedata = shapedata * R;
Modified: trunk/octave-forge/main/mechanics/inst/@rigidbody/rigidbody.m
===================================================================
--- trunk/octave-forge/main/mechanics/inst/@rigidbody/rigidbody.m 2011-10-17 20:26:22 UTC (rev 8769)
+++ trunk/octave-forge/main/mechanics/inst/@rigidbody/rigidbody.m 2011-10-18 06:38:04 UTC (rev 8770)
@@ -17,7 +17,7 @@
## @deftypefn {Function File} {@var{obj} =} rigidbody ()
## Create object @var{obj} of the rigid_body class.
##
-## If no input argument is provided the object is empty.
+## If no input argument is provided the object is filled with default values.
##
## @end deftypefn
@@ -68,40 +68,30 @@
## Shape given
if pargiven(1)
shapedata = varargin{idx(1)+1};
- if isnumeric (shapedata)
+ if !iscell(shapedata)
+ shapedata = polygon2shape(shapedata);
+ else
+ error('rigidBody:IvalidArgument','Unrecognized shape data.');
+ end
- ## Fill Shape struct
- ## Shape is always describd with axis with lower moment in positive
- ## x-axis direction and centered in [0 0].
- baricenter = center_mass_poly2d(shapedata,1);
+ ## Fill Shape struct
+ ## Shape is always describd with axis with lower moment in positive
+ ## x-axis direction and centered in [0 0].
+ baricenter = masscenter(shapedata);
- shapedata = bsxfun(@minus, shapedata, baricenter);
+ shapedata = shapetranslate(shapedata, -baricenter);
- [PA l] = principalaxes(shapedata);
- ## Put 1st axis positive in x
- if PA(1,1) < 0 ;
- PA = -PA;
- end
- rigidbody.Shape.PrincipalAxes = PA;
- rigidbody.Shape.AreaMoments = l;
-
- rigidbody.Shape.Data = varargin{idx(1)+1} * PA.';
-
- if size(shapedata,1) <= 4
- rigidbody.InertiaMoment = ...
- inertia_moment_poly2d (shapedata, rigidbody.Mass);
- else
- rigidbody.InertiaMoment = ...
- inertia_moment_ncpoly2d (shapedata, rigidbody.Mass);
- end
-
- elseif iscell (shapedata)
- #TODO
- error('rigidBody:Devel','Polycurve shape, not yet implemented.');
- else
- error('rigidBody:IvalidArgument','Unrecognized shape data.');
+ [PA l] = principalaxes(shapedata);
+ ## Put 1st axis positive in x
+ if PA(1,1) < 0 ;
+ PA = -PA;
end
+ rigidbody.Shape.PrincipalAxes = PA;
+ rigidbody.Shape.AreaMoments = l;
+ rigidbody.Shape.Data = varargin{idx(1)+1} * PA.';
+
+ rigidbody.InertiaMoment = inertiamoment (shapedata, rigidbody.Mass);
end
## Mass given
Added: trunk/octave-forge/main/mechanics/inst/deprecated/area_poly2d.m
===================================================================
--- trunk/octave-forge/main/mechanics/inst/deprecated/area_poly2d.m (rev 0)
+++ trunk/octave-forge/main/mechanics/inst/deprecated/area_poly2d.m 2011-10-18 06:38:04 UTC (rev 8770)
@@ -0,0 +1,52 @@
+%% Copyright (c) 2011 Juan Pablo Carbajal <car...@if...>
+%%
+%% This program is free software: you can redistribute it and/or modify
+%% it under the terms of the GNU General Public License as published by
+%% the Free Software Foundation, either version 3 of the License, or
+%% any later version.
+%%
+%% This program is distributed in the hope that it will be useful,
+%% but WITHOUT ANY WARRANTY; without even the implied warranty of
+%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+%% GNU General Public License for more details.
+%%
+%% You should have received a copy of the GNU General Public License
+%% along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} @var{A} = area_poly2d (@var{p})
+%% Calculates the Area of a 2D polygon.
+%%
+%% The polygon is described in @var{p}, where each row is a different vertex.
+%% The algorithm was adapted from P. Bourke web page
+%% @uref{http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/}
+%%
+%% @seealso{inertia_moment_poly2d, center_mass_poly2d}
+%% @end deftypefn
+
+function A = area_poly2d(poly)
+
+ N = size(poly,1);
+ nxt = [2:N 1];
+ px = poly(:,1);
+ px_nxt = poly(nxt,1);
+ py = poly(:,2);
+ py_nxt = poly(nxt,2);
+
+ A = sum(px.*py_nxt - px_nxt.*py)/2;
+
+end
+
+%!demo
+%! % A parametrized arbitrary triangle and its area
+%!
+%! triangle = @(a,b,h) [0 0; b 0; a h];
+%! h = linspace(0.1,1,10);
+%! b = pi;
+%! for i=1:length(h);
+%! P = triangle(0.1,b,h(i));
+%! area(i) = area_poly2d(P);
+%! end
+%!
+%! % The area of the triangle is b*h/2
+%! plot(h,area,'o',h,b*h/2);
Added: trunk/octave-forge/main/mechanics/inst/deprecated/center_mass_poly2d.m
===================================================================
--- trunk/octave-forge/main/mechanics/inst/deprecated/center_mass_poly2d.m (rev 0)
+++ trunk/octave-forge/main/mechanics/inst/deprecated/center_mass_poly2d.m 2011-10-18 06:38:04 UTC (rev 8770)
@@ -0,0 +1,60 @@
+%% Copyright (c) 2011 Juan Pablo Carbajal <car...@if...>
+%%
+%% This program is free software: you can redistribute it and/or modify
+%% it under the terms of the GNU General Public License as published by
+%% the Free Software Foundation, either version 3 of the License, or
+%% any later version.
+%%
+%% This program is distributed in the hope that it will be useful,
+%% but WITHOUT ANY WARRANTY; without even the implied warranty of
+%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+%% GNU General Public License for more details.
+%%
+%% You should have received a copy of the GNU General Public License
+%% along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} @var{cm} = center_mass_poly2d (@var{p})
+%% Calculates the center of mass a 2D polygon.
+%%
+%% The polygon is described in @var{p}, where each row is a different vertex.
+%% The algorithm was adapted from P. Bourke web page
+%% @url{http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/}
+%%
+%% @seealso{inertia_moment_poly2d, area_poly2d}
+%% @end deftypefn
+
+function cm = center_mass_poly2d(poly)
+
+ N = size(poly,1);
+ nxt = [2:N 1];
+ px = poly(:,1);
+ px_nxt = poly(nxt,1);
+ py = poly(:,2);
+ py_nxt = poly(nxt,2);
+
+ cm = zeros(1,2);
+ cr_prod = (px.*py_nxt - px_nxt.*py);
+
+ % Area
+ A = sum(cr_prod)/2;
+
+ % Center of mass
+ cm(1) = sum( (px+px_nxt) .* cr_prod );
+ cm(2) = sum( (py+py_nxt) .* cr_prod );
+ cm = cm/(6*A);
+
+end
+
+%!demo
+%! % The center of mass of this two triangles is the same
+%! % since both describe the same figure.
+%!
+%! P = [0 0; 1 0; 0 1];
+%! P2=[0 0; 0.1 0; 0.2 0; 0.25 0; 1 0; 0 1];
+%! [center_mass_poly2d(P) center_mass_poly2d(P2)],
+%!
+%! % The centroid does not give the right answer
+%!
+%! [mean(P) mean(P2)],
+
Added: trunk/octave-forge/main/mechanics/inst/deprecated/inertia_moment_ncpoly2d.m
===================================================================
--- trunk/octave-forge/main/mechanics/inst/deprecated/inertia_moment_ncpoly2d.m (rev 0)
+++ trunk/octave-forge/main/mechanics/inst/deprecated/inertia_moment_ncpoly2d.m 2011-10-18 06:38:04 UTC (rev 8770)
@@ -0,0 +1,118 @@
+%% Copyright (c) 2011 Juan Pablo Carbajal <car...@if...>
+%%
+%% This program is free software: you can redistribute it and/or modify
+%% it under the terms of the GNU General Public License as published by
+%% the Free Software Foundation, either version 3 of the License, or
+%% any later version.
+%%
+%% This program is distributed in the hope that it will be useful,
+%% but WITHOUT ANY WARRANTY; without even the implied warranty of
+%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+%% GNU General Public License for more details.
+%%
+%% You should have received a copy of the GNU General Public License
+%% along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} @var{I} = inertia_moment_ncpoly2d (@var{p}, @var{m})
+%% Calculates the moment of inertia of a simple 2D polygon.
+%%
+%% The polygon is described in @var{p}, where each row is a different vertex.
+%% @var{m} is the total mass of the polygon, assumed uniformly distributed.
+%% The polygon is triangulated using delaunay algorithm and then the
+%% Superposition Principle and the Parallel axis theorem are applied to each
+%% triangle.
+%%
+%% The position of the vertices is assumed to be given from the center of mass
+%% of the polygon.
+%% To change a general polygon to this description you can use:
+%% @code{P = P - repmat(center_mass_poly2d(P),size(P,1))}.
+%%
+%% @seealso{inertia_moment_poly2d, center_mass_poly2d}
+%% @end deftypefn
+
+function I = inertia_moment_ncpoly2d (poly, M)
+ %% Get total area
+ A = area_poly2d (poly);
+
+ %% triangulate
+ T = delaunay (poly(:,1), poly(:,2));
+ nT = size(T,1);
+
+% debug
+% triplot(T,poly(:,1),poly(:,2),'color',[0.8 0.8 0.8])
+% hold on
+% drawPolygon(poly,'k');
+% hold off;
+
+ I = 0;
+ for it = 1:nT
+ P = poly (T(it,:), :);
+ %% get centers of mass
+ cm = center_mass_poly2d (P);
+
+ % Check if triangle CoM is inside polygon
+ if ~inpolygon (cm(1), cm(2), poly(:,1), poly(:,2))
+ continue
+ end
+
+ %% get the mass as fraction of total area
+ a = area_poly2d (P);
+ if a < 0
+ aux = P(1,:);
+ P(1,:) = P(2,:);
+ P(2,:) = aux;
+ a = -a;
+ end
+ m = M*a/A;
+
+% debug
+% patch(P(:,1),P(:,2),'facecolor','b','edgecolor','none');
+% pause
+
+
+
+ %% get moment of inertia
+ mo = inertia_moment_poly2d (P, m, cm);
+
+ %% Assemble: Superposition + parallel axis
+ I += mo + m * sumsq (cm);
+ end
+
+end
+
+%!demo
+%! % C shape polygon
+%! poly = [0 0; 1 0; 1 0.25; 0.25 0.25; 0.25 0.75; 1 0.75; 1 1; 0 1];
+%!
+%! % Take to center of mass
+%! poly = poly - repmat(center_mass_poly2d(poly),size(poly,1),1);
+%! A = area_poly2d(poly);
+%!
+%! I = inertia_moment_ncpoly2d(poly,1)
+%!
+%! % It should give (breaking C in rectangles)
+%! r1 = [poly(1:3,:); poly(1,1) poly(3,2)]; a1 = area_poly2d(r1); m1 = abs(a1/A);
+%! c1 = center_mass_poly2d(r1); I1 = inertia_moment_poly2d(r1,m1,c1);
+%! r2 = [r1(4,:); poly(4:5,:); poly(1,1) poly(5,2)]; a2 = area_poly2d(r2); m2 = abs(a2/A);
+%! c2 = center_mass_poly2d(r2); I2 = inertia_moment_poly2d(r2,m2,c2);
+%! r3 = [poly(5:8,:); r2(4,:)]; a3 = area_poly2d(r3); m3 = abs(a3/A);
+%! c3 = center_mass_poly2d(r3); I3 = inertia_moment_poly2d(r3,m3,c3);
+%!
+%! I1 + m1*sumsq(c1) + I2 + m2*sumsq(c2) + I3 + m3*sumsq(c3)
+
+%!test
+%! poly = [0 0; 1 0; 1 0.25; 0.25 0.25; 0.25 0.75; 1 0.75; 1 1; 0 1];
+%! poly = poly - repmat(center_mass_poly2d(poly),size(poly,1),1);
+%! A = area_poly2d(poly);
+%! I = inertia_moment_ncpoly2d(poly,1)
+%! % It should give (breaking C in rectangles)
+%! r1 = [poly(1:3,:); poly(1,1) poly(3,2)]; a1 = area_poly2d(r1); m1 = abs(a1/A);
+%! c1 = center_mass_poly2d(r1); I1 = inertia_moment_poly2d(r1,m1,c1);
+%! r2 = [r1(4,:); poly(4:5,:); poly(1,1) poly(5,2)]; a2 = area_poly2d(r2); m2 = abs(a2/A);
+%! c2 = center_mass_poly2d(r2); I2 = inertia_moment_poly2d(r2,m2,c2);
+%! r3 = [poly(5:8,:); r2(4,:)]; a3 = area_poly2d(r3); m3 = abs(a3/A);
+%! c3 = center_mass_poly2d(r3); I3 = inertia_moment_poly2d(r3,m3,c3);
+%! I_shouldbe = I1 + m1*sumsq(c1) + I2 + m2*sumsq(c2) + I3 + m3*sumsq(c3);
+%! assert(I,I_shouldbe)
+
Added: trunk/octave-forge/main/mechanics/inst/deprecated/inertia_moment_poly2d.m
===================================================================
--- trunk/octave-forge/main/mechanics/inst/deprecated/inertia_moment_poly2d.m (rev 0)
+++ trunk/octave-forge/main/mechanics/inst/deprecated/inertia_moment_poly2d.m 2011-10-18 06:38:04 UTC (rev 8770)
@@ -0,0 +1,69 @@
+%% Copyright (c) 2011 Juan Pablo Carbajal <car...@if...>
+%%
+%% This program is free software: you can redistribute it and/or modify
+%% it under the terms of the GNU General Public License as published by
+%% the Free Software Foundation, either version 3 of the License, or
+%% any later version.
+%%
+%% This program is distributed in the hope that it will be useful,
+%% but WITHOUT ANY WARRANTY; without even the implied warranty of
+%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+%% GNU General Public License for more details.
+%%
+%% You should have received a copy of the GNU General Public License
+%% along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} @var{I} = inertia_moment_poly2d (@var{p}, @var{m}, @var{offset})
+%% Calculates the moment of inertia of a 2D star-shaped polygon.
+%%
+%% The polygon is described in @var{p}, where each row is a different vertex.
+%% @var{m} is the total mass of the polygon, assumed uniformly distributed.
+%% The optional argument @var{offset} is an origin translation vector. All vertex
+%% are transformed to the reference frame with origin at @var{offset}.
+%%
+%% This expression assumes that the polygon is star-shaped. The position of the
+%% vertices is assumed to be given from the center of mass of the polygon.
+%% To change a general polygon to this description you can use:
+%% @code{P = P - repmat(center_mass_poly2d(P),size(P,1))}.
+%% or call the function using the offset:
+%% @code{inertia_moment_poly2d (@var{p}, @var{m}, center_mass_poly2d(@var{p}))}
+%%
+%% @seealso{inertia_moment_ncpoly2d, center_mass_poly2d}
+%% @end deftypefn
+
+function I = inertia_moment_poly2d(poly,mass,offset=[])
+
+ numVerts = size(poly,1);
+
+ if !isempty (offset)
+ V = bsxfun(@minus, poly, offset);
+ Vnext = bsxfun(@minus, poly([2:numVerts 1],:),offset);
+ else
+ V = poly;
+ Vnext = poly([2:numVerts 1],:);
+ end
+
+ %% Area of the parallelograms
+ A = sqrt(sumsq(cross([Vnext zeros(numVerts,1)],[V zeros(numVerts,1)],2),2));
+ %% Distance between points
+ B = dot(V,V,2) + dot(V,Vnext,2) + dot(Vnext,Vnext,2);
+
+ C = sum(A.*B);
+
+ I = mass*C/(6*sum(A));
+end
+
+%!demo
+%!
+%! % The same triangle respect to one of its vertices described with two polygons
+%! P = [0 0; 1 0; 0 1];
+%! P2=[0 0; 0.1 0; 0.2 0; 0.25 0; 1 0; 0 1];
+%!
+%! % Now described from the center of mass of the triangle
+%! Pc = P - repmat(center_mass_poly2d(P), 3, 1);
+%! inertia_moment_poly2d(Pc,1)
+%!
+%! Pc = P2 -repmat(center_mass_poly2d(P2), size(P2,1), 1);
+%! inertia_moment_poly2d(Pc,1)
+
Added: trunk/octave-forge/main/mechanics/inst/deprecated/second_moment_poly2d.m
===================================================================
--- trunk/octave-forge/main/mechanics/inst/deprecated/second_moment_poly2d.m (rev 0)
+++ trunk/octave-forge/main/mechanics/inst/deprecated/second_moment_poly2d.m 2011-10-18 06:38:04 UTC (rev 8770)
@@ -0,0 +1,70 @@
+%% Copyright (c) 2011 Juan Pablo Carbajal <car...@if...>
+%%
+%% This program is free software: you can redistribute it and/or modify
+%% it under the terms of the GNU General Public License as published by
+%% the Free Software Foundation, either version 3 of the License, or
+%% any later version.
+%%
+%% This program is distributed in the hope that it will be useful,
+%% but WITHOUT ANY WARRANTY; without even the implied warranty of
+%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+%% GNU General Public License for more details.
+%%
+%% You should have received a copy of the GNU General Public License
+%% along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} @var{J} = second_moment_shape2d (@var{p})
+%% @deftypefnx{Function File} @var{J} = second_moment_shape2d (@var{p}),@var{type}, @var{matrix})
+%% Calculates the second moment of area of a 2D shape.
+%%
+%% The polygon is described in @var{p}, where each row is a different vertex as seen
+%% from the center of mass of the shape (see @code{center_mass_shape2d}).
+%%
+%% The output @var{J} contains Ix, Ixy and Iy, in that order. It is assumed to be
+%% oriented counter clockwise, otherwise multiply output by -1.
+%%
+%% If @var{matrix} is @code{true} then @var{J} is the symmetric second moment of area
+%% matrix, instead of the vector @code{vech(J)}.
+%% @var{type} indicates how is the shape described. So far the only case handled
+%% is when @var{shape} is a 2D polygon, where each row of @var{shape} is a vertex.
+%%
+%% @seealso{inertia_moment_shape2d, center_mass_shape2d}
+%% @end deftypefn
+
+function J = second_moment_poly2d(shape, typep='polygon', matrix=false)
+
+ if strcmpi (typep, 'polygon')
+
+ [N dim] = size (shape);
+ nxt = [2:N 1];
+ px = shape(:,1);
+ px_nxt = shape(nxt,1);
+ py = shape(:,2);
+ py_nxt = shape(nxt,2);
+
+ cr_prod = px.*py_nxt - px_nxt.*py;
+
+ J = zeros (3, 1);
+ J(1) = sum ( (py.^2 + py.*py_nxt + py_nxt.^2) .* cr_prod); % Jx
+ J(3) = sum ( (px.^2 + px.*px_nxt + px_nxt.^2) .* cr_prod); % Jy
+ J(2) = 0.5 * sum ( (px.*py_nxt + 2*(px.*py + px_nxt.*py_nxt) + px_nxt.*py) .* cr_prod);
+ J = J/12;
+
+ elseif strcmpi (typep, 'polyline')
+
+ error('second_moment_poly2d:Devel','Polyline, not implemented yet.');
+ #TODO
+ elseif strcmpi (typep, 'cbezier')
+
+ error('second_moment_poly2d:Devel', 'Cubic bezier, not implemented yet');
+ #TODO
+
+ end
+
+ if matrix
+ J = vech2mat (J);
+ end
+
+end
+
Added: trunk/octave-forge/main/mechanics/inst/inertiamoment.m
===================================================================
--- trunk/octave-forge/main/mechanics/inst/inertiamoment.m (rev 0)
+++ trunk/octave-forge/main/mechanics/inst/inertiamoment.m 2011-10-18 06:38:04 UTC (rev 8770)
@@ -0,0 +1,73 @@
+%% Copyright (c) 2011 Juan Pablo Carbajal <car...@if...>
+%%
+%% This program is free software; you can redistribute it and/or modify
+%% it under the terms of the GNU General Public License as published by
+%% the Free Software Foundation; either version 3 of the License, or
+%% (at your option) any later version.
+%%
+%% This program is distributed in the hope that it will be useful,
+%% but WITHOUT ANY WARRANTY; without even the implied warranty of
+%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+%% GNU General Public License for more details.
+%%
+%% You should have received a copy of the GNU General Public License
+%% along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} { @var{Iz} =} inertiamoment (@var{pp}, @var{m})
+%% @deftypefnx {Function File} { [@var{Iz}, @var{error}] =} inertiamoment (@dots{})
+%% @deftypefnx {Function File} { @dots{} =} inertiamoment (@dots{}, @var{tol})
+%% Moment of intertia of a plane shape.
+%%
+%% The shape is defined with piecewise smooth polynomials. @var{pp} is a
+%% cell where each elements is a 2-by-(poly_degree+1) matrix containing px(i,:) =
+%% pp{i}(1,:) and py(i,:) = pp{i}(2,:).
+%%
+%% @end deftypefn
+
+function Iz = inertiamoment (pp, m, tol=sqrt(eps)*[1 1])
+
+ Iz = sum(cellfun (@Iint, pp));
+ A = shapearea(pp);
+ Iz = m*Iz / A;
+
+endfunction
+
+function dI = Iint (x)
+
+ px = x(1,:);
+ py = x(2,:);
+
+ paux = conv (px, px)/3 + conv(py, py);
+ paux2 = conv (px, polyderiv (py)) ;
+ P = polyint (conv (paux, paux2));
+
+ %% Faster than quad
+ dI = diff(polyval(P,[0 1]));
+
+endfunction
+
+%!demo % non-convex bezier shape
+%! weirdhearth ={[-17.6816 -34.3989 7.8580 3.7971; ...
+%! 15.4585 -28.3820 -18.7645 9.8519]; ...
+%! [-27.7359 18.1039 -34.5718 3.7878; ...
+%! -40.7440 49.7999 -25.5011 2.2304]};
+%! Iz = inertiamoment (weirdhearth,1)
+
+%!test
+%! square = {[1 -0.5; 0 -0.5]; [0 0.5; 1 -0.5]; [-1 0.5; 0 0.5]; [0 -0.5; -1 0.5]};
+%! Iz = inertiamoment (square,1);
+%! assert (1/6, Iz, sqrt(eps));
+
+%!test
+%! circle = {[1.715729 -6.715729 0 5; ...
+%! -1.715729 -1.568542 8.284271 0]; ...
+%! [1.715729 1.568542 -8.284271 0; ...
+%! 1.715729 -6.715729 0 5]; ...
+%! [-1.715729 6.715729 0 -5; ...
+%! 1.715729 1.568542 -8.284271 0]; ...
+%! [-1.715729 -1.568542 8.284271 0; ...
+%! -1.715729 6.715729 0 -5]};
+%! Iz = inertiamoment (circle,1);
+%! assert (Iz, 0.5*5^2, 5e-3);
+
Added: trunk/octave-forge/main/mechanics/inst/masscenter.m
===================================================================
--- trunk/octave-forge/main/mechanics/inst/masscenter.m (rev 0)
+++ trunk/octave-forge/main/mechanics/inst/masscenter.m 2011-10-18 06:38:04 UTC (rev 8770)
@@ -0,0 +1,72 @@
+%% Copyright (c) 2011 Juan Pablo Carbajal <car...@if...>
+%%
+%% This program is free software; you can redistribute it and/or modify
+%% it under the terms of the GNU General Public License as published by
+%% the Free Software Foundation; either version 3 of the License, or
+%% (at your option) any later version.
+%%
+%% This program is distributed in the hope that it will be useful,
+%% but WITHOUT ANY WARRANTY; without even the implied warranty of
+%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+%% GNU General Public License for more details.
+%%
+%% You should have received a copy of the GNU General Public License
+%% along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} { @var{cm} =} masscenter (@var{pp}, @var{mass})
+%% @deftypefnx {Function File} { [@var{cm}, @var{error}] =} masscenter (@dots{})
+%% @deftypefnx {Function File} { @dots{} =} masscenter (@dots{}, @var{tol})
+%% Center of mass of a plane shape.
+%%
+%% The shape is defined with piecewise smooth polynomials. @var{pp} is a
+%% cell where each elements is a 2-by-(poly_degree+1) matrix containing px(i,:) =
+%% pp{i}(1,:) and py(i,:) = pp{i}(2,:).
+%%
+%% @end deftypefn
+
+function cm = masscenter (shape)
+
+ cm = sum( cell2mat ( cellfun (@CMint, shape, 'UniformOutput', false)));
+ A = shapearea(shape);
+ cm = cm / A;
+
+endfunction
+
+function dcm = CMint (x)
+
+ px = x(1,:);
+ py = x(2,:);
+ Px = polyint (conv(conv (px , px)/2 , polyderiv (py)));
+ Py = polyint (conv(-conv (px , px)/2 , polyderiv (px)));
+
+ dcm = zeros (1,2);
+ dcm(1) = diff(polyval(Px,[0 1]));
+ dcm(2) = diff(polyval(Py,[0 1]));
+
+endfunction
+
+%!demo % non-convex bezier shape
+%! weirdhearth ={[-17.6816 -34.3989 7.8580 3.7971; ...
+%! 15.4585 -28.3820 -18.7645 9.8519]; ...
+%! [-27.7359 18.1039 -34.5718 3.7878; ...
+%! -40.7440 49.7999 -25.5011 2.2304]};
+%! CoM = masscenter (weirdhearth)
+
+%!test
+%! square = {[1 -0.5; 0 -0.5]; [0 0.5; 1 -0.5]; [-1 0.5; 0 0.5]; [0 -0.5; -1 0.5]};
+%! CoM = masscenter (square);
+%! assert (CoM, [0 0], sqrt(eps));
+
+%!test
+%! circle = {[1.715729 -6.715729 0 5; ...
+%! -1.715729 -1.568542 8.284271 0]; ...
+%! [1.715729 1.568542 -8.284271 0; ...
+%! 1.715729 -6.715729 0 5]; ...
+%! [-1.715729 6.715729 0 -5; ...
+%! 1.715729 1.568542 -8.284271 0]; ...
+%! [-1.715729 -1.568542 8.284271 0; ...
+%! -1.715729 6.715729 0 -5]};
+%! CoM = masscenter (circle);
+%! assert (CoM , [0 0], 5e-3);
+
Modified: trunk/octave-forge/main/mechanics/inst/principalaxes.m
===================================================================
--- trunk/octave-forge/main/mechanics/inst/principalaxes.m 2011-10-17 20:26:22 UTC (rev 8769)
+++ trunk/octave-forge/main/mechanics/inst/principalaxes.m 2011-10-18 06:38:04 UTC (rev 8770)
@@ -14,22 +14,19 @@
## along with this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{axes} @var{l} @var{moments}] =} principalaxes (@var{shape}, type)
+## @deftypefn {Function File} {[@var{axes} @var{l} @var{moments}] =} principalaxes (@var{shape})
## Calculates the principal axes of a shape.
##
## Returns a matrix @var{axes} where each row corresponds to one of the principal
## axes of the shape. @{l} is the second moment of area around to the correspoding axis.
## @var{axes} is order from lower to higher @var{l}.
##
-## @var{type} indicates how is the shape described. So far the only case handled
-## is when @var{shape} is a 2D polygon, where each row of @var{shape} is a vertex.
-##
## @seealso{second_moment_poly2d}
## @end deftypefn
-function [PA l Jm] = principalaxes (shape, typep='polygon')
+function [PA l Jm] = principalaxes (shape)
- Jm = second_moment_poly2d (shape, typep);
+ Jm = shapemoment (shape);
Jsq = Jm(2)^2;
if Jsq > eps;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|