\\ Copyright 2015, 2016, 2017 Kevin Ryde
\\
\\ This file 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, or (at your option) any later
\\ version.
\\
\\ This file 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 file. See the file COPYING. If not, see
\\ .
\\ http://user42.tuxfamily.org/pari-geom/index.html
\\ Usage: read("geom-inertia.gp");
\\ geom_inertia_of_point([3,4]) \\ == 17/2
\\
\\
\\ This is some functions for calculating mass moment of inertia, mostly for
\\ 2-dimensional shapes in the XY plane, but some 3D. Designed for GP 2.7 up.
\\
\\ Ix = geom_inertia_Ix(M) elements of an inertia matrix
\\ Iy = geom_inertia_Iy(M)
\\ Iz = geom_inertia_Iz(M)
\\ Ixy = geom_inertia_Ixy(M)
\\ Ixz = geom_inertia_Ixz(M)
\\ Iyz = geom_inertia_Iyz(M)
\\ Ip = geom_inertia_I_at_pdir(M,p) inertia about given axis direction
\\ Ia = geom_inertia_I_at_xyangle(M,p)
\\
\\ M = geom_inertia_negate_x(M) transform axes etc of matrix
\\ M = geom_inertia_negate_y(M)
\\ M = geom_inertia_negate_z(M)
\\ M = geom_inertia_transpose_xy(M)
\\ M = geom_inertia_rotate_xydir(M,p)
\\
\\ M = geom_inertia_of_point(p) construct for various shapes
\\ M = geom_inertia_of_points(v)
\\ M = geom_inertia_of_segment(p1,p2)
\\ M = geom_inertia_of_path(v)
\\ M = geom_inertia_of_rectangle(w,h)
\\ M = geom_inertia_of_rectangle_m1(w,h)
\\ M = geom_inertia_of_triangle(p1,p2)
\\ M = geom_inertia_of_triangle_m1(p1,p2)
\\ M = geom_inertia_of_polygon(v)
\\ M = geom_inertia_of_intformal_rows(ymin,ymax,xmin,xmax,density)
\\ M = geom_inertia_of_circle(r)
\\ M = geom_inertia_of_circle_m1(r)
\\ M = geom_inertia_of_ellipse(r)
\\ M = geom_inertia_of_ellipse_m1(r)
\\
\\ p = geom_inertia_principal_xydir(M)
\\ a = geom_inertia_principal_xyangle(M)
\\ p = geom_inertia_principal_2xydir(M)
\\
\\ Mass moment of inertia for rotation about a given axis is
\\
\\ Ia = sum(m * r^2)
\\
\\ where the sum is over points of mass m at perpendicular distance r from
\\ the axis. Similarly integral for a continuous shape. Inertia is ratio
\\ of torque to angular acceleration for a rigid body rotating about this
\\ axis, similar to the way ordinary mass is the ratio of force to linear
\\ acceleration.
\\
\\ An inertia tensor has inertia for rotation about three axes X, Y, Z, and
\\ also cross products which are used to give inertia about other axis lines.
\\ An inertia tensor is represented here as a matrix
\\
\\ [ Ix, -Ixy, -Ixz ]
\\ M = [ -Iyx, Iy, -Iyz ]
\\ [ -Izx, -Izy, Iz ]
\\
\\ The matrix elements are sums (or integrals) over each point p=[x,y,z] in
\\ a shape
\\
\\ Ix = sum(y^2 + z^2) Ixy = Iyx = sum(x*y)
\\ Iy = sum(x^2 + z^2) Ixz = Iyx = sum(x*z)
\\ Iz = sum(x^2 + y^2) Iyz = Iyx = sum(y*z)
\\
\\ For example, Ix is sum of squared perpendicular distances from p to the X
\\ axis. This is the inertia when the shape rotates around the X axis.
\\
\\ The code here is mostly for plane figures which is z=0 for all points.
\\ This simplifies to
\\
\\ [ Ix, -Ixy, 0 ] Ix = sum(y^2) Ixy = Iyx = sum(x*y)
\\ M = [ -Iyx, Iy, 0 ] Iy = sum(x^2)
\\ [ 0, 0, Iz ] Iz = Ix + Iy
\\
\\ Notice that Ix is sum over y^2. Remember this is for rotation around the
\\ X axis, so y distances are wanted. Vice versa for Iy.
\\
\\ The Ixy etc cross-product terms are sum(x*y) etc. They are negated in
\\ the inertia matrix. Taking "Ixy" to mean plain sum(x*y), rather than the
\\ negative in the matrix, is intended to be easiest to remember if forming
\\ a cross-product explicitly. These products are positive for a shape
\\ entirely in 1st or 3rd quadrants, or negative if there's more sum in the
\\ 2nd or 4th quadrant.
\\
\\ The matrix form allows axes to be rotated by multiplying a matrix of
\\ rotation. Suppose R is a 3x3 matrix of rotation which rotates
\\ coordinates so point x,y,z becomes x',y',z'.
\\
\\ [ x' ] [ x ] R = rotation matrix
\\ [ y' ] = R * [ y ] matdet(R)=1
\\ [ y' ] [ z ] R^-1 = mattranspose(R)
\\
\\ then inertia matrix M' for the shape in these new coordinates is
\\
\\ M' = R^-1 * M * R
\\
\\ This is a rotation of the axes by R. If the axes are to be held fixed
\\ and the shape rotated then use the inverse of R (so rotate the axes the
\\ other way). For rotating the shape it can help to think where the axes
\\ should end up in the rotated shape.
\\
\\
\\ Shape Construction:
\\
\\ geom_inertia_of_polygon() is the most general plane shape function. The
\\ triangle and rectangle functions come out the same but take arguments in
\\ different style.
\\
\\ Inertia of a complicated shape can be constructed by adding inertia of
\\ separate pieces. For holes in a shape the inertia of the holes can be
\\ subtracted.
\\
\\ The main functions are for uniform distribution of mass, with total mass
\\ equal to length of a line or area of a shape. The "m1" functions are
\\ instead mass 1 uniformly distributed in the shape. For simple shapes the
\\ m1 forms are the basic formulas and are multiplied by area etc for
\\ density of mass=area etc.
\\
\\ Inertia is often calculated for rotation about axes passing through the
\\ centroid (centre of gravity) of a shape. Pieces with axes through the
\\ centroid can be shifted by the parallel axis theorem. To shift the
\\ origin to point p simply add
\\
\\ m * geom_inertia_of_point(p)
\\
\\ where m is the mass of the shape. This can arise when constructing a
\\ shape out of other parts.
\\
\\ If axes are not through the centroid then a shape can be shifted by first
\\ subtracting m*geom_inertia_of_point(centroid) to take axes back to the
\\ centroid, then add to the desired new location. Nothing in the inertia
\\ matrix records centroid, so this must be known separately.
\\
\\ Some plane shapes with a non-uniform distribution of mass can be handled
\\ by geom_inertia_of_intformal_rows(). It's intended for use with
\\ polynomials to describe the shape and mass distribution. A polynomial
\\ approximation to something more complex might be put through. The
\\ integrals it does are not particularly complicated, but the function is a
\\ handy way to have them all at once. (Something similar over a line
\\ parameterized by position t along its length or a 3D shape would be
\\ possible.)
\\
\\ Many shape construction functions don't enquire much into their x,y point
\\ values, so can take points which have further variables such as 'w or 'h
\\ for some overall size etc. It's best to do this with the vector point
\\ forms, since not certain polynomials with complex or quad coefficients
\\ quite work properly yet.
\\
\\ \\ triangle with vertices at a,b etc
\\ geom_inertia_of_polygon(['a,'b], ['c,'d], ['e,'f])
\\
\\
\\ Principal Axes:
\\
\\ Any 3x3 real symmetric matrix (such as inertia) can be rotated to
\\ diagonalize it. The eigenvectors of the matrix are "principal axes" of
\\ the shape which is where cross products Ixy=Ixz=Iyz=0. The physical
\\ significance is that rotation about an axis with cross products 0 is
\\ perfectly balanced, exerting no torque on the mount points as it rotates.
\\
\\ If a shape has some symmetry then there might be many principal axes.
\\ For example a sphere is the same in all directions so every direction is
\\ principal.
\\
\\ When inertia is measured over directions p the result is an ellipsoid.
\\ This can be seen in the formula of geom_inertia_I_at_pdir(). The
\\ semi-axes of the ellipsoid are the principal axes of the shape.
\\
\\ Various shapes can have the same inertia matrix. A rectagular or
\\ elliptical shape can be chosen, oriented to the principal axes and with
\\ suitable size, to reproduce any given inertia matrix. There's nothing
\\ here to do that, but it's not too hard to work back from the formulas
\\ after finding principal axes.
\\
\\ For a shape entirely in the XY plane, the Z axis through the centroid is
\\ always a principal axis, for rotation within the XY plane. Principal
\\ axes for a shape in the plane rotating up out of the plane can be found
\\ by geom_inertia_principal_xydir(). Roughly speaking these axes are where
\\ there is equal quantity of shape each side of the axis line, weighted by
\\ r^2 perpendicular distance to the axis.
\\
\\ geom_inertia_principal_xydir() is a simple 2x2 matrix eigenvector finder.
\\ It starts from geom_inertia_principal_2xydir() which is a double-angle
\\ direction. If you have quads in the inertia components then the halved
\\ xydir might not be exact. In that case 2xydir and axis direction as
\\ 1/2*atan(y/x) might be preferred.
\\
\\
\\----------------------------------------------------------------------------
\\ Version 1 - the first version.
\\-----------------------------------------------------------------------------
if(geom_point_x == 'geom_point_x, read("geom.gp")); \\ if not already loaded
geom_inertia_of_Ix_Iy_Ixy(Ix,Iy,Ixy) =
{
[ Ix, -Ixy, 0;
-Ixy, Iy, 0;
0, 0, Ix+Iy ];
}
{addhelp(geom_inertia_of_Ix_Iy_Ixy,
"M = geom_inertia_of_Ix_Iy_Ixy(Ix,Iy,Ixy)
Ix = sum(y^2), Iy=sum(x^2), Ixy=sum(x*y) are components of inertia
for a shape in the XY plane.
Return an inertia matrix for these components.
Ix and Iy are inertia rotating about the X or Y axes respectively. With
coordinate z=0 for a shape in the XY plane they are Ix=sum(y^2) and
Iy=sum(x^2).")};
geom_inertia_Ix(M) = M[1,1];
geom_inertia_Iy(M) = M[2,2];
geom_inertia_Iz(M) = M[3,3];
geom_inertia_Ixy(M) = - M[1,2];
geom_inertia_Ixz(M) = - M[1,3];
geom_inertia_Iyz(M) = - M[2,3];
\\
{addhelp(geom_inertia_Ix,
"Ix = geom_inertia_Ix(M)
M is a 3x3 inertia matrix.
Return its inertia for rotating around the X axis.")};
{addhelp(geom_inertia_Iy,
"Iy = geom_inertia_Iy(M)
M is a 3x3 inertia matrix.
Return its inertia for rotating around the Y axis.")};
{addhelp(geom_inertia_Iz,
"Iz = geom_inertia_Iz(M)
M is a 3x3 inertia matrix.
Return its inertia for rotating around the Z axis.")};
{addhelp(geom_inertia_Ixy,
"Ixy = geom_inertia_Ixy(M)
M is a 3x3 inertia matrix.
Return its XY cross-product of inertia Ixy=sum(x*y).
(If X and Y are principal axes then Ixy is 0.)")};
\\------------------------------------------------------------------------------
\\ Transformations
geom_inertia_negate_x(M) = \\ negate x so mirror across YZ plane
{
[ M[1,1], -M[1,2], -M[1,3];
-M[2,1], M[2,2], M[2,3];
-M[3,1], M[3,2], M[3,3] ];
}
{addhelp(geom_inertia_negate_x,
"M = geom_inertia_negate_x(M)
M is a 3x3 inertia matrix.
Return a new inertia matrix for the shape with X coordinates negated,
ie. mirror image across the YZ plane.")};
geom_inertia_negate_y(M) = \\ negate y so mirror across XZ plane
{
[ M[1,1], -M[1,2], M[1,3]; \\ x*x y*x z*x
-M[2,1], M[2,2], -M[2,3]; \\ x*y y*y
M[3,1], -M[3,2], M[3,3] ]; \\ x*z z*z
}
{addhelp(geom_inertia_negate_y,
"M = geom_inertia_negate_y(M)
M is a 3x3 inertia matrix.
Return a new inertia matrix for the shape with Y coordinates negated,
ie. mirror image across the XZ plane.")};
geom_inertia_negate_z(M) = \\ negate z so mirror across XY plane
{
[ M[1,1], M[1,2], -M[1,3];
M[2,1], M[2,2], -M[2,3];
-M[3,1], -M[3,2], M[3,3] ];
}
{addhelp(geom_inertia_negate_z,
"M = geom_inertia_negate_z(M)
M is a 3x3 inertia matrix.
Return a new inertia matrix for the shape with all Z coordinates negated,
ie. mirror image across the XY plane.
For a shape entirely in the XY plane this is no change at all.")};
geom_inertia_transpose_xy(M) = \\ swap 1<->2
{
[ M[2,2], M[2,1], M[2,3];
M[1,2], M[1,1], M[1,3];
M[3,2], M[3,1], M[3,3] ];
}
{addhelp(geom_inertia_transpose_xy,
"M = geom_inertia_transpose_xy(M)
M is a 3x3 inertia matrix.
Return a new inertia matrix for the shape with all X,Y coordinates
transposed (swapped), ie. mirror image across the plane defined by the
leading diagonal X=Y and by the Z axis.")};
geom_inertia_rotate_xydir(M,p) =
{
my(x = geom_point_x(p),
y = geom_point_y(p),
xx = x^2,
yy = y^2,
xy = x*y,
n = xx + yy,
Ix = geom_inertia_Ix(M),
Iy = geom_inertia_Iy(M),
Ixy = geom_inertia_Ixy(M));
geom_inertia_of_Ix_Iy_Ixy(( Ix*xx - 2*Ixy*xy + Iy*yy )/n,
( Iy*xx + 2*Ixy*xy + Ix*yy )/n,
( (Ix-Iy)*xy + Ixy*(xx-yy) )/n);
}
{addhelp(geom_inertia_rotate_xydir,
"M = geom_inertia_rotate_xydir(M,p)
M is a 3x3 inertia matrix of a shape entirely in the XY plane.
p is a point in the XY plane defining a direction.
Return a new inertia matrix which is for coordinates rotated so the new X
axis is in the direction of point p. Only the direction of p is used, its
magnitude is ignored. For example p=1+I rotates the XY axes by +45 degrees.
If you want to hold the axes fixed and rotate the shape then apply the
opposite rotation (negate the Y part of p). Or think of where the X axis
should be after the shape rotates, and give that direction for p.")};
\\-----------------------------------------------------------------------------
\\ Points
geom_inertia_of_point(p) =
{
my(x=geom_point_x(p),
y=geom_point_y(p),
z=geom_point_z(p),
xx=x^2, yy=y^2, zz=z^2,
Ixy=x*y, Ixz=x*z, Iyz=y*z);
[ yy+zz, -Ixy, -Ixz;
-Ixy, xx+zz, -Iyz;
-Ixz, -Iyz, xx+yy ];
}
{addhelp(geom_inertia_of_point,
"M = geom_inertia_of_point(p)
Return an inertia matrix for a unit mass at point p (anywhere 3-D).
Multiply the result to change the mass if desired.")};
geom_inertia_of_points(v) = sum(i=1,#v, geom_inertia_of_point(v[i]));
{
addhelp(geom_inertia_of_points,
"M = geom_inertia_of_points(v)
v is a vector of points. Return an inertia matrix for all the points.
Each point is taken to have mass 1. Multiply the result to give them all
another mass if desired.
This is simply the sum of geom_inertia_of_point(z) for each point in v.")};
\\-----------------------------------------------------------------------------
\\ Lines
geom_inertia_of_segment_m1(p1,p2) = \
geom_inertia_of_point(p1-p2)/12 + geom_inertia_of_point((p1+p2)/2);
{
addhelp(geom_inertia_of_segment_m1,
"M = geom_inertia_of_segment_m1(p1,p2)
Return an inertia matrix for a line segment between points p1 and p2
(anywhere 3-D).
The segment is taken to have mass 1 uniformly distributed along its length.
Multiply the result to change that mass if desired.
If p1==p2 then the result is mass 1 at that location, the same as
geom_inertia_of_point().
The inertia is 1/12*geom_inertia_of_point(p1-p2) for line centred at the
origin, plus geom_inertia_of_point((p1+p2)/2) to move to the actual midpoint
by the parallel axis theorem. If p1==p2 then the first term is 0, giving
point mass at p1=p2.")};
geom_inertia_of_segment(p1,p2) = \
geom_inertia_of_segment_m1(p1,p2) * geom_length_of_segment(p1,p2);
{
addhelp(geom_inertia_of_segment,
"M = geom_inertia_of_segment(p1,p2)
Return an inertia matrix for a line segment between points p1 and p2
(anywhere 3-D).
The segment is taken to have mass equal to its length and uniformly
distributed along the length. Multiply the result to change this mass if
desired.
If p1==p2 then the length is 0 so mass 0 and inertia all 0s.")};
geom_inertia_of_path(v) = \
sum(i=1,#v-1, geom_inertia_of_segment(v[i], v[i+1]));
{
addhelp(geom_inertia_of_path,
"M = geom_inertia_of_path(v)
v is vector of points which are a path of points (anywhere 3-D).
Return an inertia matrix of the whole path. The path is taken to have mass
equal to its length, uniformly distributed along its length. Multiply the
result to change the mass if desired.
This is a sum of segment inertias
geom_inertia_of_segment(v[1],v[2])
+ ...
+ geom_inertia_of_segment(v[#v-1],v[#v])
If v has just 1 point or is empty then it is taken as zero-length and
inertia matrix all 0s.
No attention is paid to whether the segments overlap. If they do then those
sections will have masses doubled (etc).")};
geom_inertia_of_path_m1(v) =
{
my(l=geom_length_of_path(v));
if(l,
geom_inertia_of_path(v) / l,
geom_inertia_of_point(v[1]));
}
{addhelp(geom_inertia_of_path_m1,
"M = geom_inertia_of_path_m1(v)
v is vector of points which are a path of points (anywhere 3-D).
Return an inertia matrix for the whole whole path. The path is taken to
have total mass 1, uniformly distributed along its length. Multiply the
result to change this mass if desired.
v must contain at least 1 point, since otherwise the location of the mass is
unspecified. If the path length is zero, either due to the path being a
single point, or the same point repeated, then the inertia is mass 1 at that
location, the same as geom_inertia_of_point().")};
\\-----------------------------------------------------------------------------
\\ Rectangles
geom_inertia_of_rectangle_m1(w,h) =
{
geom_inertia_of_Ix_Iy_Ixy(1/12*h^2,
1/12*w^2,
0);
}
{addhelp(geom_inertia_of_rectangle_m1,
"M = geom_inertia_of_rectangle_m1(w,h)
Return the inertia matrix for a rectangle width w and height h in the XY
plane centred on the origin. Should have w,h >= 0. The rectangle is taken
to have mass 1, uniformly distributed.
If h=0 then the return is the same as a geom_inertia_of_segment_m1() from
-w/2 to +w/2. Similarly if w=0 a vertical line +/-h/2. If both w=h=0 then
it is inertia all 0s (the same as geom_inertia_of_point([0,0])).")};
geom_inertia_of_rectangle(w,h) = \
geom_inertia_of_rectangle_m1(w,h) * w*h;
{
addhelp(geom_inertia_of_rectangle,
"M = geom_inertia_of_rectangle(w,h)
Return the inertia matrix of a rectangle width w and height h in the XY
plane centred at the origin. Should have w,h >= 0. The rectangle is taken
to have mass equal to its area w*h, uniformly distributed.")};
\\-----------------------------------------------------------------------------
\\ Triangle
\\ + +
\\ | |
\\ | z2 =c,d | z2 z1
\\ | * | *------*
\\ | / \ | / __-
\\ | / __-* z1 =a,b | / __--
\\ |/ - |/ -
\\ *---------+--- *---------+---
\\
\\ c.f. area = 1/2 * (a*d - b*c);
\\
geom_inertia_of_triangle_m1(p1,p2,p3=0) =
{
(geom_inertia_of_point(p1+p2)
+ geom_inertia_of_point(p2+p3)
+ geom_inertia_of_point(p3+p1) ) / 12;
}
{addhelp(geom_inertia_of_triangle_m1,
"M = geom_inertia_of_triangle_m1(p1,p2,p3=0)
Return an inertia matrix for a triangle with vertices at p1, p2, p3
(anywhere 3-D). p3 defaults to the origin 0. The vertices can be given in
any order.
The triangle is taken to have mass 1 uniformly distributed over its area.
Multiply the result to change this mass if desired.
This inertia is equivalent to mass 1/3 at the midpoint of each triangle
side.
If p1==p2==p3 then the result is mass 1 at that location, the same as
geom_inertia_of_point().
If p1,p2,p3 are colinear then note that the result is not the same as a line
segment between the ends. Effectively the sloping sides of the triangle are
a non-uniform distribution of mass along the length, which is maintained
down to a colinear limit.")};
geom_inertia_of_triangle(p1,p2,p3=0) = \
geom_inertia_of_triangle_m1(p1,p2,p3) * geom_area_of_triangle(p1,p2,p3);
{
addhelp(geom_inertia_of_triangle,
"M = geom_inertia_of_triangle(p1,p2,p3=0)
Return an inertia matrix for a triangle with vertices at p1, p2, p3 in the
XY plane. p3 defaults to the origin 0.
The triangle is taken to have mass equal to its area, uniformly distributed.
Multiply the result to change that mass if desired.
The vertices should go anti-clockwise around (the usual mathematical
direction). If they go clockwise then the area is negative, in the manner
of geom_area_of_triangle(), so the mass is negative.")};
\\-----------------------------------------------------------------------------
\\ Polygon
geom_inertia_of_polygon(v) = \
sum(i=1,#v, geom_inertia_of_triangle(v[i], if(i==#v,v[1],v[i+1])));
{
addhelp(geom_inertia_of_polygon,
"M = geom_inertia_of_polygon(v)
v is a polygon in the form of a vector of points in the XY plane.
Return an inertia matrix for area inside the polygon.
The polygon is taken to have mass equal to its area, uniformly distributed.
The points should go anti-clockwise around (the usual mathematical
direction). If they go clockwise then the area is negative as per
geom_area_of_polygon() and so the mass is negative.")};
\\-----------------------------------------------------------------------------
\\ Circle
geom_inertia_of_circle_m1(r) = \
my(Ir=1/4*r^2); geom_inertia_of_Ix_Iy_Ixy(Ir,Ir,0);
{
addhelp(geom_inertia_of_circle_m1,
"M = geom_inertia_of_circle_m1(r)
Return the inertia of a circle of radius r centred at the origin in the XY
plane. The circle is taken to have mass 1, uniformly distributed. If r=0
then the inertia is all 0s, the same as for a point mass at the origin.
The inertia is Ix = Iy = 1/4*r^2, with cross-product Ixy=0 since
symmetric.")};
geom_inertia_of_circle(r) = \
my(Ir=Pi/4*r^4); geom_inertia_of_Ix_Iy_Ixy(Ir,Ir,0);
{
addhelp(geom_inertia_of_circle,
"M = geom_inertia_of_circle(r)
Return the inertia of a circle radius r at the origin in the XY plane.
The circle is taken to have mass equal to its area Pi*r^2, uniformly
distributed.")};
\\-----------------------------------------------------------------------------
\\ Ellipse
geom_inertia_of_ellipse_m1(a,b) = \
geom_inertia_of_Ix_Iy_Ixy(1/4*b^2, 1/4*a^2, 0);
{
addhelp(geom_inertia_of_ellipse_m1,
"M = geom_inertia_of_ellipse_m1(a,b)
Return the inertia of a ellipse in the XY plane centred at the origin and
aligned to the XY axes. a is the radius along the X axis. b is the radius
along the Y axis. The ellipse is taken to have mass 1, uniformly
distributed.
The inertia is Ix=1/4*b^2 and Iy=1/4*a^2, with cross-product Ixy=0 since
symmetric in mirror image across the X or Y axes (every point x*y has a
corresponding -x*y, for net sum Ixy=0).
If a=b=0 then it is the same as a point mass 1 at the origin, which is
inertia 0.
If one of a,b=0 then note that the result is not the same as a line between
the ends. Effectively the elliptical shape is a non-uniform distribution of
mass along the length, which is maintained down to a line limit.")};
geom_inertia_of_ellipse(a,b) = \
Pi*a*b * geom_inertia_of_ellipse_m1(a,b);
{
addhelp(geom_inertia_of_ellipse,
"M = geom_inertia_of_ellipse(r)
Return the inertia of a ellipse centred at the origin in the XY plane.
a is the radius along the X axis, b is the radius along the Y axis. The
ellipse is taken to have mass equal to its area Pi*a*b, uniformly
distributed.
If a=0 or b=0 then the area is 0 and the inertia is 0.")};
\\ ellipse area pi*a*b a=X radius, b=Y radius
\\ ellipse Ix = 1/4 * pi*a*b^3
\\ ellipse Iy = 1/4 * pi*a^3*b
\\ ellipsoid volume 4/3*pi*a*b*c (?)
\\ ellipsoid Ix = 1/5 * mass * (b^2+c^2)
\\ ellipsoid Iy = 1/5 * mass * (a^2+c^2)
\\ ellipsoid Iz = 1/5 * mass * (a^2+b^2)
\\-----------------------------------------------------------------------------
\\ intformals
geom_inertia_of_intformal_rows(ymin,ymax, xmin,xmax, density=1) =
{
my(double_intformal =
(poly)->
poly = intformal(poly*density, 'x);
poly = subst(poly,'x,xmax) - subst(poly,'x,xmin);
poly = intformal(poly,'y);
subst(poly,'y,ymax) - subst(poly,'y,ymin));
geom_inertia_of_Ix_Iy_Ixy(double_intformal('y*'y), \\ Ix
double_intformal('x*'x), \\ Iy
double_intformal('x*'y)); \\ Ixy
}
{addhelp(geom_inertia_of_intformal_rows,
"M = geom_inertia_of_intformal_rows(ymin,ymax, xmin,xmax)
Return the inertia matrix for shape in the XY plane defined by an integral
over rows.
A given row goes from x=xmin to x=xmax. xmin,xmax can be any polynomial in
the Y coordinate variable 'y, so allowing sides like lines or parabolas.
density is a polynomial in variables 'x,'y representing mass density as a
function of location x,y. The default density=1 gives mass proportional to
area and uniformly distributed across the area.
Rows are integrated from y=ymin to y=ymax.
y=ymax x=xmax
Ix = integral integral density*y^2 dx dy
y=ymin x=xmin
similarly Iy = integral of density*x^2
and Ixy = integral of density*x*y
The integrals are made with intformal() so all of xmin,xmax, ymin,ymax and
density can be polynomials in further variables too, for example some 'w or
'h as overall width and height. This gives a parameterized final formula
and is a typical usage. A final formula might then be massaged manually to
make an inertia calculating function (or just subst() into it).
See geom_area_of_intformal_rows() to calculate area. For simple shapes
(like rectangles or triangles) the area may be a divisor of the inertia
formula.")};
\\-----------------------------------------------------------------------------
geom_inertia_I_at_pdir(M,p) =
{
my(x = geom_point_x(p),
y = geom_point_y(p),
z = geom_point_z(p));
( geom_inertia_Ix(M)*x^2
+ geom_inertia_Iy(M)*y^2
+ geom_inertia_Iz(M)*z^2
- 2*( x*y*geom_inertia_Ixy(M)
+ x*z*geom_inertia_Ixz(M)
+ y*z*geom_inertia_Iyz(M)) )
/ geom_norm(p);
}
{addhelp(geom_inertia_I_at_pdir,
"Ip = geom_inertia_I_at_pdir(M,p)
M is an inertia matrix. Return inertia Ip of the shape rotating about an
axis in the direction of point p (anywhere 3-D).
The formula is
( Ix*x^2 + Iy*y^2 + Iz*z^2 - 2*(x*y*Ixy + x*z*Ixz + y*z*Iyz) )
/ (x^2 + y^2 + z^2)
If p is on the X axis (so y=z=0) then the return is just Ix, similarly Iy
or Iz. The magnitude of p is ignored, just its direction is used.
As a function of direction p, the result is an ellipsoid oriented to the
principal axes of M. For p in the XY plane, so z=0, it reduces to an
ellipse
(Ix*x^2 - 2*x*y*Ixy + Iy*y^2 ) / (x^2 + y^2).")};
geom_inertia_I_at_xyangle(M,alpha) =
{
my(c=cos(alpha), s=sin(alpha));
geom_inertia_Ix(M)*c^2 - 2*geom_inertia_Ixy(M)*c*s + geom_inertia_Iy(M)*s^2;
}
{
addhelp(geom_inertia_I_at_xyangle,
"Ia = geom_inertia_I_at_xyangle(M, alpha)
M is an inertia matrix. Return the inertia Ia of the shape rotating about
an axis in the XY plane at angle alpha radians anti-clockwise around from
the X axis. This is
Ia = Ix*cos^2 - 2*Ixy*cos*sin + Iy*sin^2")};
\\
\\ Ix + Iy
\\ I(alpha) = ------- + D*cos(2*alpha - delta)
\\ 2
\\
\\ D = abs(d) delta=arg(d)
\\
\\ Ix - Iy
\\ d = ------- - Ixy * I
\\ 2
\\
\\ minimum
\\ 2*alpha-delta = Pi
\\ 2*alpha = delta + Pi
\\ maximum with cos(2*alpha-delta) = 1
\\ 2*alpha-delta = 0
\\ 2*alpha = delta
\\
geom_inertia_principal_2xydir(M) = \
[ geom_inertia_Ix(M) - geom_inertia_Iy(M), -2*geom_inertia_Ixy(M) ];
{
addhelp(geom_inertia_principal_2xydir,
"point = geom_inertia_principal_2xydir(M)
M is an inertia matrix for a shape in the XY plane.
Return a point p which is at an angle from the X axis twice the principal
axis with maximum inertia.
The magnitude of p is unspecified, just its direction should be used.
-p is a point at angle twice the principal axis with minimum inertia.")};
geom_inertia_principal_xydir(M) = \
geom_xydir_half(geom_inertia_principal_2xydir(M));
{
addhelp(geom_inertia_principal_xydir,
"point = geom_inertia_principal_xydir(M)
M is an inertia matrix for a shape in the XY plane.
Return a point p which is on the principal axis with maximum inertia.
The magnitude of p is unspecified, just its direction should be used.
This is geom_xydir_half() of geom_inertia_principal_2xydir(). If there are
quads in M is a quad then it's possible the halving won't be exact. In that
case the suggestion would be use geom_inertia_principal_2xydir() and
principal axis is at 1/2*atan(y/x) of that.
The prinipal axis with minimum inertia is at +90 degrees from p. That is
[-y,x], eg. using geom_point_rotate_xydir(p,[0,1]).")};
geom_inertia_principal_xyangle(M) =
{
my(p=geom_inertia_principal_xydir(M));
arg(1.0*geom_point_x(p) + I*1.0*geom_point_y(p));
}
{addhelp(geom_inertia_principal_xyangle,
"angle = geom_inertia_principal_xyangle(M)
M is an inertia matrix for a shape in the XY plane.
Return the angle in radians of the principal axis where inertia is the
maximum. The angle is in the range (-Pi/2, +Pi/2].
The prinipal axis where inertia is the minimum is at right angle to the
maximum, so +/-Pi/2 from a.")};
\\------------------------------------------------------------------------------
\\ LocalWords: Ryde geom gp XY Iy Iz Ixy xy xydir Ixz Iyx Iyz Izx Izy YZ
\\ LocalWords: matdet mattranspose addhelp ie XZ intformals intformal dx
\\ LocalWords: ymin ymax xmin xmax poly subst dy parameterized yy zz cos
\\ LocalWords: xbase abs arg min arctan pdir RFRACs diagonalizes colinear