Updates

Author: melund

InstantaneousAxisOfRotation.any

@@ -57,40 +57,58 @@ The object will IAOR object will plot the axis, and provide the following output
    AnyRefFrame& Body2,
    SHOW_LINEAR_VELOCITY = "On",
    AXIS_LOCATION = "MIDPOINT",
-   SHOW_SHORTEST_DIST = "On"
-){
+   SHOW_SHORTEST_DIST = "On")
+{
 
+  AnyKinLinear Linear ={
+    Ref=0;
+    AnyRefFrame  &Body1= .Body1;
+    AnyRefFrame  &Body2= .Body2;
+  };
+  
+  AnyKinRotational Rotational ={
+    Type = RotVector;
+    AngVelOnOff = On;
+    AnyRefFrame  &Body1= .Body1;
+    AnyRefFrame  &Body2= .Body2;
+  };  
+   
+  /// Size of the drawn spheres
   #var AnyVar SphereSize=0.02;
+  
+  /// Length of the axis drawn
   #var AnyVar AxisLength=2;
+  
+  /// Thickness of the axis drawn
   #var AnyVar AxisThickness = 0.01;
 
 
-     /// Rotational Velocity of Body2 with repect to Body1
-     AnyVec3 omega = Rotational.Vel;
-     
-     /// Linear Velociy of Body2 with respect to Body1
-     AnyVec3 r_dot = Linear.Vel; 
-     
-     /// Postion of Body2 with respect to Body1
-     AnyVec3 r = Linear.Pos;
-     
-     /// The rotation magnitude 
-     AnyVar RotationMagnitude = iffun(eqfun(vnorm(omega),0),1e-15, vnorm(omega));
-     
-     /// The rotation axis direction 
-     AnyVec3 RotationAxisDirection = omega/RotationMagnitude;
-     
-     /// The point on the rotation axis closest to Body2
-     AnyVec3 rc_2 = r + cross(omega, r_dot)/(RotationMagnitude^2);
-     
-     /// The point on the rotation axis closest to Body1
-     AnyVec3 rc_1 = rc_2-rc_2*RotationAxisDirection'*RotationAxisDirection;
-     
-     /// The motion pitch or the ratio of linear motion to angular motion
-     AnyVar MotionPitch = (omega*r_dot')/(RotationMagnitude^2);
-     
-     /// Linear Velocity of points on the rotation axis
-     AnyVec3 rc_dot = MotionPitch*omega;
+  /// Rotational Velocity of Body2 with repect to Body1
+  AnyVec3 omega = Rotational.Vel;
+  
+  /// Linear Velociy of Body2 with respect to Body1
+  AnyVec3 r_dot = Linear.Vel; 
+  
+  /// Postion of Body2 with respect to Body1
+  AnyVec3 r = Linear.Pos;
+  
+  /// The rotation magnitude 
+  AnyVar RotationMagnitude = iffun(eqfun(vnorm(omega),0),1e-15, vnorm(omega));
+  
+  /// The rotation axis direction 
+  AnyVec3 RotationAxisDirection = omega/RotationMagnitude;
+  
+  /// The point on the rotation axis closest to Body2
+  AnyVec3 rc_2 = r + cross(omega, r_dot)/(RotationMagnitude^2);
+  
+  /// The point on the rotation axis closest to Body1
+  AnyVec3 rc_1 = rc_2-rc_2*RotationAxisDirection'*RotationAxisDirection;
+  
+  /// The motion pitch or the ratio of linear motion to angular motion
+  AnyVar MotionPitch = (omega*r_dot')/(RotationMagnitude^2);
+  
+  /// Linear Velocity of points on the rotation axis
+  AnyVec3 rc_dot = MotionPitch*omega;
 
      AnyFolder Visualizations = 
      {
@@ -169,17 +187,6 @@ The object will IAOR object will plot the axis, and provide the following output
      };
      #endif
    };
-   AnyKinLinear Linear ={
-     Ref=0;
-     AnyRefFrame  &Body1= .Body1;
-     AnyRefFrame  &Body2= .Body2;
-   };
-   
-   AnyKinRotational Rotational ={
-     Type = RotVector;
-     AngVelOnOff = On;
-     AnyRefFrame  &Body1= .Body1;
-     AnyRefFrame  &Body2= .Body2;
-   };   
+ 
 
 };

Readme.md

@@ -1,44 +1,115 @@
 # Class for calculating Instantaneous Axis of Rotation
 
-The `InstantaneousAxisOfRotation` class template will calculate the instantenous axis of rotation between two bodies. 
+The instantaneous axis of rotation between two bodies is a useful concept in
+biomechanics. In this post, we will dig into how to calculate the instantaneous
+axis of rotation and show off an AnyScript `class_template` to do this
+calculation and plot the axis between any two reference frames in the AnyBody
+Modeling System.
 
-Instantenous axis of rotation are the points fixed to a body undergoing planar movement that has zero velocity at a particular instant of time. At this instant, the velocity vectors of the trajectories of other points in the body generate a circular field around this point which is identical to what is generated by a pure rotation.
+The are several papers that discusses the
+importance of the instantaneous axis of rotation in biomechanics. See for
+example XX and XX.  But first, let us look briefly at the math behind.
 
-Planar movement of a body is often described using a plane figure moving in a two-dimensional plane. The instant centre is the point in the moving plane around which all other points are rotating at a specific instant of time.
+## 3D Rotations
 
-The continuous movement of a plane has an instant centre for every value of the time parameter. This generates a curve called the moving centrode. The points in the fixed plane corresponding to these instant centres form the fixed centrode.
+We can view any displacement of a body in a three-dimensional space as rotation
+around a single axis and a translation along that axis. This gives rise to the
+idea of a screw motion in 3D along what is also called the helical axis. If we
+spit the movement up into infinitesimally small movements each of these will
+have an rotation axis and a linear velocity along that axis. This is the
+instantaneous axis of rotation, and as the name indicate will change direction
+and location as the object moves. 
 
-------------------------------------------------------------------------------
+If we know the angular and linear velocity of any point on a rigid body, we can
+calculate the properties of the screw motion and the instantaneous axis of
+rotation.
 
-Think of rotation as a manifestation of a change in coordinate direction. Affix a coordinate system at _any_ point on a rigid body and it is going to change direction _at the same rate_ regardless of the point location (even at the rotation axis). This is why angular velocity is shared among the entire rigid body.
+From a mathematical point of view the rate of change of a 3x3 rotation matrix
+$R$ is:
 
-Mathematically, the rate of change of a 3×3 rotation matrix $R$ is $$ \frac{{\rm d}}{{\rm d}t}R = \vec{\omega} \times R$$ where $\vec{\omega}$ is the angular rotation vector, and $\times$ the vector cross product.
+$$ \frac{{\mathrm d}}{{\mathrm d}t}R = \vec{\omega} \times R$$
 
-For the second part of your question a general rotation _can_ be decomposed into three rotations about three predefined axes (like precession, nutation and spin), or by a single rotation about an arbitrary axis. In either case, the rotation is going to be about a single axis and you need three parameters to specify it exactly. Either three scalar angle speeds about three predefined axes, or two parameters for the rotation axis, and one for the rotation speed magnitude.
+where $\vec{\omega}$ is called the angular velocity vector. This angular
+velocity is quantity shared by all points on the rigid body (including those on
+the rotation axis). Using $\vec{\omega}$ we can write the linear velocity
+$\vec{v}$ of any point at some distance $\vec{r} from the instantaneous axis of
+rotation as:
 
-Note also that the general motion of a rigid body is a rotation about some axis in 3D, coupled with a parallel translation along the same axis. The rotation and translation defines a screw in 3D and given the linear and angular velocities of _any_ point on the rigid body, the properties of the motion can be derived.
+$$\vec{v} = \vec{\omega} \times \vec{r}$$ (1)
 
-For example, said rigid body has rotational velocity $\vec{\omega}$ and linear velocity $\vec{v}_A$ at some point _A_ with position $\vec{r}_A$ at some instant. The following properties are defined:
+Similarly we can also calculate backwards and find the intantanous axis of
+rotation at some distance $-\vec{r}$ from any point if we know the linear
+velocity of the point and the angular velocity:
 
- 1. The rotation axis direction is $$\vec{e} = \frac{\vec{\omega}}{|\vec{\omega}|}$$
- 2. The rotation magnitude is $$\omega = |\vec{\omega}|$$
- 3. The point _C_ on the rotation line closest to the origin is $$\vec{r}_C = \vec{r}_A + \frac{\vec{\omega}\times\vec{v}_A}{\omega^2}$$
- 4. The motion pitch (ratio of linear motion to angular motion) is $$ h = \frac{\vec{\omega} \cdot \vec{v}_A}{\omega^2}$$
- 5. The linear velocity of point _C_ is $$\vec{v}_C = h\,\vec{\omega}$$
+This is given by:
 
-In summary, a single rotation speed, and a linear velocity vector at _some_ point is enough to describe the instantaneous motion of a rigid body both geometrically (rotation line in 3D, pitch and magnitude) and analytically (rigid body transformation of velocity to any other point).
+$$ 
+ -\vec{r}= \frac{\vec{\omega} \times \vec{v} }{ \vec{\omega}\cdot\vec{\omega} } 
+$$
 
+Here is short proof. Skip it if you like: 
 
----------------------------------------------------------
+>Start with $\vec{\omega} \times \vec{v}$, insert (1) and use the  tripple product to expand it:
+>$$ 
+>\vec{\omega} \times \vec{v} = \vec{\omega} \times (\vec{\omega} \times \vec{r})=\vec{\omega}(\vec{\omega}\cdot\vec{r})-\vec{r}(\vec{\omega}\cdot\vec>{\omega})
+>$$
+>
+>The shortest vector $\vec{r}$ from the rotation axis to any point is by defintion perpendicular to angular velocity $\vec{\omega}$, hence:
+>
+>$$ 
+>\require{cancel}\vec{\omega} \times \vec{v} =\vec{\omega}(\cancel{\vec{\omega}\cdot\vec{r}})-\vec{r}(\vec{\omega}\cdot\vec{\omega})=-r(\vec{\omega}>\cdot\vec{\omega})
+>$$
 
 
 
+Thus if know the position of a point $\vec{r}_P$ we can find closest point
+$\vec{r}_C$ on the rotation axis as:
 
+ $$\vec{r}_C = \vec{r}_P - \vec{r} = \vec{r}_P + \frac{\vec{\omega}\times\vec{v}_P}{\vec{\omega}\cdot\vec{\omega}}$$
 
 
 
+## Properties
 
+If for example we have a rigid body with angular velocity $\vec{\omega}$ and
+some point $P$ with position $\vec{r}_P$ and velocity $\vec{v}_p$ then we can
+define the different properties:
 
+1. The magintude of agular rotation:
+$$ \omega = \|\vec{\omega}\|$$
+
+2. The direction of the intantanous axis of rotation:
+
+$$ \vec{v}_{IOAR} = \frac{\vec{\omega}}{\omega}$$
+
+3. The point C on the intantanous axis of rotation:
+
+ $$\vec{r}_C = \vec{r}_P + \frac{\vec{\omega}\times\vec{v}_P}{\vec{\omega}\cdot\vec{\omega}}$$
+
+4. Ratio of angular to linear angular velocity:
+
+ $$ h = \frac{\vec{\omega} \cdot \vec{v}_P}{\vec{\omega}\cdot\vec{\omega}}$$
+
+5. The linear velocity at point $C$ along the screw axis is then given by
+
+ $$\vec{v}_C = h\,\vec{\omega}$$
+
+
+All we need to know is the rotation velocity of a body and the velocity of any point to find the instantaneous axis of rotation. 
+
+## AnyScript implementation
+
+In AnyScript we can easily find the angular velocity between two reference frames using the class `AnyKinRotational` and the linear velocity using the class `AnyKinLinear`. 
+
+We have created a custom class_template that makes it easy to add all this code to a model and display the instantaneous axis of rotation. Here is a short example on how to use the class. 
+
+
+The class template is called `InstantaneousAxisOfRotation`.  does exactly this, and then calculate the properties listed above. 
+
+
+
+
+------------------
 
 
 
@@ -47,7 +118,7 @@ In summary, a single rotation speed, and a linear velocity vector at _some_ poin
 
 
 EXAMPLE:
-```
+```c++
 #include "<ANYBODY_PATH_TOOLBOX>/Rotations/InstantaneousAxisOfRotation.any"
 
 Main = {