(c) 1993, 2009  M.Lampton
                             STELLAR SOFTWARE

     The orientation of an object in three dimensions can be
     described an a variety of ways.  One such description as
     follows: take an arbitrary point P that is not the origin, and
     write its coordinates (x,y,z) in a frame of reference fixed in
     the object.  Also, write its coordinates (X,Y,Z) in a lab frame
     of coordinates having the same center.  The matrix M that
     converts (x,y,z) into (X,Y,Z) is uniquely defined by the
     orientation of the object's coordinate frame.    Explicitly:

             X       | M11  M12  M13 |   x
             Y   =   | M21  M22  M23 |   y
             Z       | M31  M32  M33 |   z

     This description is a set of nine numbers, namely the matrix
     elements. They are not independent: the sums of the squares
     along any column or row is exactly 1.0 because the length of
     the (x,y,z) vector and the (X,Y,Z) vector have to come out
     equal for any object vector.  A matrix that satisfies these six
     constraints is said to be unitary.  Nine numbers and six
     constraints means that there are just three degrees of freedom
     for the three dimensional orientation.

     Orientations can also be described using angles.  For any
     rotation, there exists some direction in space (two parameters)
     and one angle around that axis that gives the rotation.  In
     Cartesian coordinates it is usually more convenient to break
     this one eigenaxis rotation into three separate consecutive
     rotations about the original or partly rotated axes.

     How many such descriptions are there?   The first motion can
     be taken about the lab X, Y, or Z axes -- three possibilities.
     Whichever is chosen, the second rotation cannot be around that
     axis again, but now there are two new axes to choose from, plus
     two unused of the original three axes.  If the second angle is
     one of the original three, it creates three new axes for a
     total of eight, seven of which are possible third axes, giving
     42 possible descriptions.  If the second angle is one of the
     two new axes produced by the first angle, it creates two new
     axes for a total of seven, six of which are possible third
     angle axes, giving 36 more descriptions.  In all there are 78.


     In BEAM products, and in much of the optics industry, we use
     the rotation-sequence described by successive rotations about
     the X, y', and z" axes.   This is just one way of doing
     business, but it has the advantage that each coordinate appears
     once so that if you have a simple motion about X, Y, or Z,
     there is a one-angle description of it, and it is about the
     obvious axis.  Another advantage of (X,y',z") over say (X,Y,Z)
     is that this is the way that actual mechanical goniometers
     work: each goniometer cradle carries the base of the next
     goniometer, so that the first motion is about a fixed lab
     coordinate direction -- the second motion is about an axis that
     depends on the first motion, and the third motion axis depends
     on both the preceding motions.  In all, 18 of the rotation
     sequences are of the form (A, B', C") that can be created with
     actual goniometers.  Of these, (X, y', z") is the easiest to

     Define these three angles in a positive right-handed sense:

             Tilt "t" -- about lab +X axis
             Pitch "p"  -- about tilted +y axis
             Roll "r"  -- about tilted and pitched +z axis.

     Abbreviate:     ct = cos(t)
                     st = sin(t)
                     cp = cos(p)
                     sp = sin(p)
                     cr = cos(r)
                     sr = sin(r)

     Then, any vector (xrot, yrot, zrot) seen in the frame of the
     rotated coordinates can be converted to its lab frame
     representation by the linear matrix operator M:

         Xlab  = |     cp*cr           -sr*cp         sp   |  xrot
         Ylab  = | st*sp*cr+ct*sr   ct*cr-sr*st*sp  -cp*st |  yrot
         Zlab  = | -sp*ct*cr+st*sr  sp*sr*ct+st*cr   ct*cp |  zrot

     This matrix M is unitary.

     For the case of a flat surface whose orientation is defined by
     the tilt-pitch-roll coordinates, its normal is the unit +z axis
     in its local (rotated) frame.  In lab coordinates this +z axis,
     and hence the normal, is just the rightmost column of the
     matrix M:

                     Xlab = sin(p)
                     Ylab = -cos(p)*sin(t)
                     Zlab = cos(p)*cos(t)

     The inverse of any unitary matrix M is its transpose.  The
     conversion of a lab frame vector into that vector expressed in
     the rotated coordinates is done with the transpose of M, called

         xrot  = | cp*cr   st*sp*cr+ct*sr  -sp*ct*cr+st*sr |  Xlab
         yrot  = | -sr*cp  ct*cr-sr*st*sp   sp*sr*ct+st*cr |  Ylab
         zrot  = |   sp        -cp*st            ct*cp     |  Zlab


     Another of the many three-angle descriptions is the (X,Y,Z)
     sequence:  the first rotation is identical to tilt (right
     handed about the lab X axis), the second rotation is about the
     lab Y axis even though the optic is no longer aligned with the
     lab frame, and the third rotation is about the lab Z axis,

     Call these angles ax, ay, and az, and abbreviate

             cx = cos(ax)
             sx = sin(ax)
             cy = cos(ay)
             sy = sin(ay)
             cz = cos(az)
             sz = sin(az)

     Then, any point (xrot yrot zrot) in the frame of the rotated
     coordinates can be converted to its lab frame representation by
     the linear operator  M(Z) * M(Y) * M(X) because successive
     motions in the lab frame multiply onto the left hand side (the
     lab side) of the building matrix.  This product is:

         Xlab  = | cy*cz  -cx*sz+sx*sy*cz   cx*sy*cz+sx*sz |  xrot
         Ylab  = | cy*sz   cx*cz+sx*sy*sz  -sx*cz+cx*sy*sz |  yrot
         Zlab  = | -sy        sx*cy              cx*cy     |  zrot

     There is a set of angles (ax, ay, az) that describes the same
     orientation as the set of angles (tilt, pitch, roll).  Of
     course the angles are different, but the orientations are the
     same.  This fact gives us a way to convert angles in one scheme
     to or from angles in another.  We set up the rotation matrix
     for each scheme.  All nine of the matrix elements are equal
     because the orientations are the same.  So, solve either way
     for the set of angles.   This plan converts any rotation
     description to any other.

     Here's an example connecting (ax, ay, az) to (tilt, pitch,
     roll).  The third element is sin(pitch) in the tilt-pitch-roll
     system. But at the same time the third element is also given by
     cos(ax)*sin(ay)*cos(az)+sin(ax)*sin(az).  So, we have a formula
     for pitch, given ax ay az:

         pitch = arcsin[ cos(ax)*sin(ay)*cos(az) + sin(ax)*sin(az) ]

     Similarly, note that the ratio of the sixth matrix element to
     the ninth matrix element is -tan(tilt).  So:

                         -sin(ax)*cos(az) + cos(ax)*sin(ay)*sin(az)
         tilt = -arctan[ ------------------------------------------ ]

     Similarly, the ratio of the second matrix element to the first
     matrix element is -tan(roll), and

                         -cos(ax)*sin(az) + sin(ax)*sin(ay)*cos(az)
         roll = -arctan[ ------------------------------------------ ]

     By shopping for the simplest expression of variables you need
     in one matrix, and using the equality of matrix elements, you
     can write down any needed conversion formula.

     This works in both directions.  For example let's find the
     angles ax ay and az that correspond to a given trio of tilt,
     pitch, and roll.   First, ax can be obtained by noting that
     tan(ax) is the ratio of the eighth to ninth matrix element:

                      sin(p)*sin(r)*cos(t) + sin(t)*cos(r)
         ax = arctan[ ------------------------------------ ]

     The term sin(ay) appears by itself as the seventh matrix
     element, so

         ay = -arcsin[ -sin(p)*cos(t)*cos(r) + sin(t)*sin(r) ].

     Finally, the tangent of az is the ratio of the fourth matrix
     element of M to its first element:

                       sin(t)*sin(p)*cos(r) + cos(t)*sin(r)
         az = arctan[  ------------------------------------ ]


     Yet another of the specifications is the (Y,X,Z) system.  The
     transformation M = R(Z) * R(X) * R(Y) works out to:

         Xlab  = | cy*cz-sx*sy*sz  -cx*sz   sy*cz+sx*cy*sz |  xrot
         Ylab  = | cy*sz+sx*sy*cz   cx*cz   sy*sz-sx*cy*cz |  yrot
         Zlab  = |   -cx*sy           sx         cx*cy     |  zrot

     Again, useful formulas going either way can be found by
     identifying corresponding parts of this matrix and the (t,p,r)
     matrix that performs the same transformation, from rotated to
     lab coordinates.


     This is Euler's original set of angles and is widely used in 
     dynamics although not much used in optics.  Adopt the angle
     names a = roll about Z, i = pitch about x', and p = roll about z".
     Adopt the abbreviations

             ca = cos(a)
             sa = sin(a)
             ci = cos(i)
             si = sin(i)
             cp = cos(p)
             sp = sin(p)

     Then the transformation M = R(Z) * R(x') * R(z") works out to:     

         Xlab  = | ca*cp-sa*ci*sp  -ca*sp-sa*ci*cp   sa*si |  xrot
         Ylab  = | sa*cp+ca*ci*sp   ca*ci*cp-sa*sp  -ca*si |  yrot
         Zlab  = |    si*sp             si*cp          ci  |  zrot