SHARE
TWEET

Mathematica Magnetic Field Line .gif

a guest Jul 11th, 2016 107 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* Numerically integrate by Gaussian quadrature. MUCH faster than NIntegrate. *)
  2. Needs["NumericalDifferentialEquationAnalysis`"];
  3. {qx,qweights}=Transpose[GaussianQuadratureWeights[10,0,1]];
  4. integrate[f_,a_,b_]=(b-a)qweights.f/@((b-a)qx+a);
  5. (* Simple helper function *)
  6. norm[v_]:=Sqrt[v.v];
  7. (* 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 *)
  8. 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];
  9. (* B is the sum of two fields *)
  10. B[x_]:=fieldFromLoop[x,1,0,0]+fieldFromLoop[x,1,2,Pi/7];
  11. (* numerically integrate/draw a field line, by setting r'[t]=B[r[t]] *)
  12. tmax=1000000;
  13. 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}];
  14. (* Draw two loops and a B field *)
  15. loop1=ParametricPlot3D[{Cos[theta],Sin[theta],0},{theta,0,2 Pi},PlotStyle->Black]/.Line[pts_,rest___]:>Tube[pts,0.1,rest];
  16. 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];
  17. bfield=ParametricPlot3D[{x[t],y[t],z[t]}/.soln,{t,0,tmax},MaxRecursion->15,SphericalRegion->True,Axes->False,Boxed->False];
  18. (* Export a .gif! *)
  19. 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]
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top