(*package for plotting vector fields*)

Needs["VectorFieldPlots`"]

(*Domain*)
n=4; (*number of subintervals in each direction*)
X={x,-2,2};dX=(X[[3]]-X[[2]])/n;AppendTo[X,dX];
Y={y,-2,2};dY=(Y[[3]]-Y[[2]])/n;AppendTo[Y,dY];
Z={z,-2,2};dZ=(Z[[3]]-Z[[2]])/n;AppendTo[Z,dZ];
r=Min[{dX,dY,dZ}];

(*Vector Field*)
F[x_,y_,z_]={x z,x y z,-(y^2)};
FIELD=VectorFieldPlot3D[F[x,y,z],X,Y,Z,VectorHeads->True]

(*Curl and magnitude of the curlthe vector field at a point*)
Curl[x_,y_,z_]={(∂_y)F[x,y,z][[3]]-(∂_z)F[x,y,z][[2]],(∂_z)F[x,y,z][[1]]-(∂_x)F[x,y,z][[3]],(∂_x)F[x,y,z][[2]]-(∂_y)F[x,y,z][[1]]};
Mag[v_]:=Module[{a},a=v;Simplify[(Sqrt[a.a])]]

(*basis for the plane normal to the curl vector*)
(*u[x_,y_,z_]:=((F[x,y,z])/(Mag[F[x,y,z]]));*)
v[x_,y_,z_,F_]:=Module[{a},a=Simplify[Cross[Curl[x,y,z],F]];(a/(Mag[a]))]
v[x_,y_,z_]:=v[x,y,z,F[x,y,z]]
u[x_,y_,z_]:=Cross[Curl[x,y,z],v[x,y,z]]

(*The range of values taken on by the magnitude of the curl vectors - this helps scale the graphics*)
CurlRange=Flatten[Table[Mag[Curl[i,j,k]],{i,-2,2,0.8},{j,-2,2,0.8},{k,-2,2,0.8}]];
CurlRange={Min[CurlRange],Max[CurlRange]};

(*Circle formed by the curl "eddy"*)
Eddy[x_,y_,z_,t_]=r Cos[t]u[x,y,z]+r Sin[t]v[x,y,z]+{x,y,z}//N;

(*paddle wheel spun by the curl "eddy"*)
Paddle[x_,y_,z_,t_]:=Graphics3D[Polygon[{{x,y,z},Eddy[x,y,z,t],{x,y,z}+ ((3 r)/(CurlRange[[2]]))Curl[x,y,z]}]]

Manipulate[Show[{ParametricPlot3D[Eddy[1,1,1,t],{t,0,2π}],Paddle[1,1,1,t],Paddle[1,1,1,t+((2π)/3)],Paddle[1,1,1,t-((2π)/3)]},PlotRange->4],{t,0,2π,0.01}]

Manipulate[Animate[Show[{Paddle[x,y,z,t],Paddle[x,y,z,t+((2π)/3)],Paddle[x,y,z,t-((2π)/3)],FIELD},BoxRatios->1],{t,0,2π,0.1}],{x,-2,2,0.8},{y,-2,2},{z,-2,2}]