Mathematica Magnetic Field Line .gif

(* Numerically integrate by Gaussian quadrature. MUCH faster than NIntegrate. *)
Needs["NumericalDifferentialEquationAnalysis`"];
{qx,qweights}=Transpose[GaussianQuadratureWeights[10,0,1]];
integrate[f_,a_,b_]=(b-a)qweights.f/@((b-a)qx+a);
(* Simple helper function *)
norm[v_]:=Sqrt[v.v];
(* give the magnetic field at position vector x, from a loop of radius r, centered at some weird rotated position controlled by z and phi *)
fieldFromLoop[x_,r_,z_,phi_]:=integrate[Function[theta,1/(4 Pi)Cross[{-r Sin[theta],r Cos[theta],0},rprime]/norm[rprime]^3/.rprime->((x-RotationMatrix[phi,{1,1,0}].{r Cos[theta],r Sin[theta],z}))],0,2Pi];
(* B is the sum of two fields *)
B[x_]:=fieldFromLoop[x,1,0,0]+fieldFromLoop[x,1,2,Pi/7];
(* numerically integrate/draw a field line, by setting r'[t]=B[r[t]] *)
tmax=1000000;
soln=NDSolve[Join[{x[0]==0,y[0]==0,z[0]==0},Thread[{x'[t],y'[t],z'[t]}==B[{x[t],y[t],z[t]}]]],{x,y,z},{t,0,tmax}];
(* Draw two loops and a B field *)
loop1=ParametricPlot3D[{Cos[theta],Sin[theta],0},{theta,0,2 Pi},PlotStyle->Black]/.Line[pts_,rest___]:>Tube[pts,0.1,rest];
loop2=ParametricPlot3D[RotationMatrix[phi,{1,1,0}].{r Cos[theta],r Sin[theta],z}/.{r->1,z->2,phi->Pi/7},{theta,0,2 Pi},Evaluated->True,PlotStyle->Black]/.Line[pts_,rest___]:>Tube[pts,0.1,rest];
bfield=ParametricPlot3D[{x[t],y[t],z[t]}/.soln,{t,0,tmax},MaxRecursion->15,SphericalRegion->True,Axes->False,Boxed->False];
(* Export a .gif! *)
Export["bField.gif",Table[Show[bfield,loop1,loop2,PlotRange->6,ViewPoint->({Cos[h],Sin[h],1}),ImageSize->{200,200}],{h,0,2 Pi,2 Pi/120}],"DisplayDurations"->1/30]