Euler Angles
Dave Typinski, 2010

Notes: 
1) Change ONLY the data in the User Selectable Constants section

Order of rotation:
first in α,
then in β,
then in γ.

R = Rγ Rβ Rα  = RZxz

where

    * The rightmost matrix represents the α rotation around the axis ' z ' of the original reference frame
    * The middle matrix represents the β rotation around an intermediate ' x ' axis which is the "line of nodes".
    * The leftmost matrix represents the γ rotation around the axis ' Z ' of the final reference frame

Clear[α,β,γ];

(* point p in original reference frame x,y,z *)
p={x,y,z};

(* rotation matrix for rotation about z in α *)
Rz={{Cos[α],Sin[α],0},
{-Sin[α],Cos[α],0},
{0,0,1}};

(* rotation matrix for rotation about X in β *)
Rx={{1,0,0},
{0,Cos[β],Sin[β]},
{0,-Sin[β],Cos[β]}};

(* rotation matrix for rotation about Z in γ *)
RZ={{Cos[γ],Sin[γ],0},
{-Sin[γ],Cos[γ],0},
{0,0,1}};

(* rotation matrix for rotation about z, then x, then Z *)
RZxz=RZ.Rx.Rz;
RZxz//MatrixForm

(* rotation about z, then x, then Z *)
prot=RZxz.p//MatrixForm

A={{a,b,c},{d,e,f},{g,h,i}};
B={x,y,z};
Print[A//MatrixForm,B//MatrixForm," = ",A.B//MatrixForm]

Clear[θ,ϕ,φ];

(* point p in original reference frame x,y,z *)
p={x,y,z};

(* rotation matrix for rotation about x LH *)
Rxl={{1,0,0},
{0,Cos[φ],-Sin[φ]},
{0,Sin[φ],Cos[φ]}};

(* rotation matrix for rotation about y LH *)
Ryl={{Cos[ϕ],0,Sin[ϕ]},
{0,1,0},
{-Sin[ϕ],0,Cos[ϕ]}};

(* rotation matrix for rotation about z LH *)
Rzl={{Cos[θ],-Sin[θ],0},
{Sin[θ],Cos[θ],0},
{0,0,1}};

(* rotation about x, then y, then z LH*)
Rzl.Ryl.Rxl//MatrixForm

(* rotation vector about x, then y, then z LH*)
Rzl.Ryl.Rxl.p//MatrixForm

Clear[θ,ϕ,φ];

(* point p in original reference frame x,y,z *)
p={x,y,z};

(* rotation matrix for rotation about x RH *)
Rxr={{1,0,0},
{0,Cos[φ],Sin[φ]},
{0,-Sin[φ],Cos[φ]}};

(* rotation matrix for rotation about y RH *)
Ryr={{Cos[ϕ],0,-Sin[ϕ]},
{0,1,0},
{Sin[ϕ],0,Cos[ϕ]}};

(* rotation matrix for rotation about z RH *)
Rzr={{Cos[θ],Sin[θ],0},
{-Sin[θ],Cos[θ],0},
{0,0,1}};

(* rotation about x, then y, then z RH*)
(* Rzr.Ryr.Rxr//MatrixForm *)

(* rotation vector about x, then y, then z RH*)
Rzr.Ryr.Rxr.p//MatrixForm

Clear[θ,ϕ,φ];
points={{0,0,3},{21.873,0,3},{21.873,0,3.656},{0,0,3.656}}
ListPointPlot3D[points,PlotRange->{{-5,25},{-10,10},{-10,10}},AxesLabel->{x,y,z}]
θ=5.*π/180;
ϕ=10.*π/180;
φ=30.*π/180;
points[[All,3]]=points[[All,3]]-3;
points
Do[points[[i]]=Rzr.Ryr.Rxr.points[[i]],{i,4}];
points
points[[All,3]]=points[[All,3]]+3;
points
ListPointPlot3D[points,PlotRange->{{-25,25},{-10,10},{-10,10}},AxesLabel->{x,y,z}]

(* EZNEC and Mathematica use right handed 3D coordinate systems 
we use a right handed rotation matrix to invert the sense of rotation opposite to
that commonly used in math and physics so that it's more user-friendly 
i.e., a + rotation in the XY plane rotates the +Y axis towards the +X axis *)