\\ Copyright 2015, 2016, 2017, 2018 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. Uses geom.gp. \\ addhelp() strings have docs of each function and see general notes in the \\ comments below. \\ \\ 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,a) \\ \\ 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_Ix_Iy_Ixy(Ix,Iy,Ixy) \\ 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_segment_m1(p1,p2) \\ M = geom_inertia_of_path(v) \\ M = geom_inertia_of_path_m1(v) \\ M = geom_inertia_of_rectangle(w,h) \\ M = geom_inertia_of_rectangle_m1(w,h) \\ M = geom_inertia_of_triangle(p1,p2,p3=0) \\ M = geom_inertia_of_triangle_m1(p1,p2,p3=0) \\ M = geom_inertia_of_polygon(v) \\ M = geom_inertia_of_intformal_rows(ymin,ymax, xmin,xmax, density=1) \\ M = geom_inertia_of_circle(r) \\ M = geom_inertia_of_circle_m1(r) \\ M = geom_inertia_of_ellipse(a,b) \\ M = geom_inertia_of_ellipse_m1(a,b) \\ \\ p = geom_inertia_principal_xydir(M) \\ angle = geom_inertia_principal_xyangle(M) \\ p = geom_inertia_principal_2xydir(M) \\ \\ \\ Inertia \\ ------- \\ \\ 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 distance r from the axis, with \\ that distance measured perpendicular to the axis. Similarly an 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. geom_inertia_of_Ix_Iy_Ixy() puts separate sum terms \\ into a matrix. \\ \\ \\ Rotation \\ -------- \\ \\ 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 = matrix of rotation \\ [ 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 \\ matrices 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, but probably less useful.) \\ \\ 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 (described in geom.gp), since not sure 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 are zero 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. \\ \\ If inertia is measured over various directions p then the resulting \\ magnitudes are 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 but 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. \\ \\ \\ Requirements \\ ------------ \\ \\ Designed for GP 2.7 up. Clean to option "strictargs". \\ \\ \\ Changes \\ ------- \\ \\ Version 1 - the first version. \\ Version 2 - doc tweaks. \\----------------------------------------------------------------------------- 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 = sqr(x), yy = sqr(y), 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=sqr(x), yy=sqr(y), zz=sqr(z), 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*sqr(h), 1/12*sqr(w), 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*sqr(r)); 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*sqr(b), 1/4*sqr(a), 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(a,b) 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, density=1) Return the inertia matrix for a 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)*sqr(x) + geom_inertia_Iy(M)*sqr(y) + geom_inertia_Iz(M)*sqr(z) - 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 Ip = ( 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 resulting Ip 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,a) = { my(c=cos(a), s=sin(a)); geom_inertia_Ix(M)*sqr(c) - 2*geom_inertia_Ixy(M)*c*s + geom_inertia_Iy(M)*sqr(s); } { addhelp(geom_inertia_I_at_xyangle, "Ia = geom_inertia_I_at_xyangle(M,a) M is an inertia matrix. Return the inertia Ia of the shape rotating about an axis in the XY plane at angle \"a\" radians anti-clockwise around from the X axis. This is Ia = Ix*cos(a)^2 - 2*Ixy*cos(a)*sin(a) + Iy*sin(a)^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, "p = 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, "p = 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 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 principal 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 principal 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 intformal dx \\ LocalWords: ymin ymax xmin xmax poly subst dy parameterized yy zz cos \\ LocalWords: xbase abs arg min arctan pdir RFRACs diagonalizes colinear \\ LocalWords: Ip Ia xyangle atan eg